{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving a two-asset HANK model in sequence space\n",
    "\n",
    "In this notebook we solve the two-asset HANK model from Auclert, Bardóczy, Rognlie, Straub (2019): \"Using the Sequence-Space Jacobian to Solve and Estimate Heterogeneous-Agent Models\". Please see the [Github repository](https://github.com/shade-econ/sequence-jacobian) for more information and code.\n",
    "\n",
    "The only new concept relative to the one-asset HANK notebook is that of a **solved block**. These are extensions of simple blocks that are useful for large macro models and that's why we introduce them here, in the context of our richest model. We will apply all other concepts and tools that are familiar from other notebooks here without further explanation. \n",
    "\n",
    "## Model description\n",
    "\n",
    "The household problem is characterized by the Bellman equation\n",
    "\n",
    "$$\n",
    "\\begin{align} \\tag{1}\n",
    "V_t(e_t, b_{t-1}, a_{t-1}) = \\max_{c_t, b_t, a_t} &\\left\\{\\frac{c_t^{1-\\sigma}}{1-\\sigma} + \\beta E_t V_{t+1}(e_{t+1}, b_t, a_t) \\right\\}\n",
    "\\\\\n",
    "c_t + a_t + b_t &= (1-\\tau_t)w_t n_t e_t + (1 + r_t^a)a_{t-1} + (1 + r_t^b)b_{t-1} - \\Psi(a_t, a_{t-1}) \n",
    "\\\\\n",
    "a_t &\\geq 0, \\quad b_t \\geq \\underline{b},\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where the adjustment cost function is specified as\n",
    "\n",
    "$$\n",
    "\\Psi(a_t, a_{t-1}) =  \\frac{\\chi_1}{\\chi_2}\\left|\\frac{a_t - (1 + r_t^a) a_{t-1}}{(1 + r_{t}^a) a_{t-1} + \\chi_0}\\right|^{\\chi_2} \\left[(1 + r_t^a) a_{t-1} + \\chi_0 \\right],\n",
    "$$\n",
    "\n",
    "with $\\chi_0, \\chi_1 > 0$ and $\\chi_2 > 1.$ For the full description of the model, including the problems of the other agents, please see appendix A.3 of the paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0 Import packages\n",
    "\n",
    "The first two are standard python packages, the rest contain code we wrote for this project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import utils\n",
    "import jacobian as jac\n",
    "from het_block import het\n",
    "from simple_block import simple\n",
    "from solved_block import solved\n",
    "import two_asset\n",
    "import nonlinear\n",
    "import determinacy as det"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Calibrate steady state\n",
    "\n",
    "We developed an efficient backward iteration function to solve the Bellman equation in (1). Although we view this as a contribution on its own, discussing the algorithm goes beyond the scope of this notebook. If you are interested in how we solve a two-asset model with convex portfolio-adjustment costs in discrete time, please see appendix B of the paper for a detailed description and `two_asset.py` for the implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ss = two_asset.hank_ss(noisy=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Define simple blocks\n",
    "\n",
    "Compare these to the equations in appendix A.3 of the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "@simple\n",
    "def dividend(Y, w, N, K, pi, mup, kappap, delta):\n",
    "    psip = mup / (mup - 1) / 2 / kappap * np.log(1 + pi) ** 2 * Y\n",
    "    I = K - (1 - delta) * K(-1)\n",
    "    div = Y - w * N - I - psip\n",
    "    return psip, I, div\n",
    "\n",
    "@simple\n",
    "def taylor(rstar, pi, phi):\n",
    "    i = rstar + phi * pi\n",
    "    return i\n",
    "\n",
    "@simple\n",
    "def fiscal(r, w, N, G, Bg):\n",
    "    tax = (r * Bg + G) / w / N\n",
    "    return tax\n",
    "\n",
    "@simple\n",
    "def finance(i, p, pi, r, div, omega, pshare):\n",
    "    rb = r - omega\n",
    "    ra = pshare * (div + p) / p(-1) + (1-pshare) * (1 + r) - 1\n",
    "    fisher = 1 + i(-1) - (1 + r) * (1 + pi)\n",
    "    return rb, ra, fisher\n",
    "\n",
    "@simple\n",
    "def wage(pi, w, N, muw, kappaw):\n",
    "    piw = (1 + pi) * w / w(-1) - 1\n",
    "    psiw = muw / (1 - muw) / 2 / kappaw * np.log(1 + piw) ** 2 * N\n",
    "    return piw, psiw\n",
    "\n",
    "@simple\n",
    "def union(piw, N, tax, w, U, kappaw, muw, vphi, frisch, beta):\n",
    "    wnkpc = kappaw * (vphi * N**(1+1/frisch) - muw*(1-tax)*w*N*U) + beta * np.log(1 + piw(+1)) - np.log(1 + piw)\n",
    "    return wnkpc\n",
    "\n",
    "@simple\n",
    "def mkt_clearing(p, A, B, Bg):\n",
    "    asset_mkt = p + Bg - B - A\n",
    "    return asset_mkt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Define solved blocks\n",
    "\n",
    "Solved blocks are mini SHADE models embedded as blocks inside larger SHADE models. Like simple blocks, solved blocks correspond to aggregate equilibrium conditions: they map sequences of aggregate inputs directly into sequences of aggregate outputs. The difference is that in the case of simple blocks, this mapping has to be analytical, while solved blocks are designed to accommodate implicit relationships that can only be evaluated numerically. \n",
    "\n",
    "Such implicit mappings between variables become more common as macro complexity increases. Solved blocks are a valuable tool to simplify the DAG of large macro models.\n",
    "\n",
    "### 3.1 Price setting \n",
    "The Phillips curve characterizes $(\\pi)$ conditional on $(Y, mc, r):$ \n",
    "\n",
    "$$\n",
    "\\log(1+\\pi_t) = \\kappa_p \\left(mc_t - \\frac{1}{\\mu_p} \\right) + \\frac{1}{1+r_{t+1}} \\frac{Y_{t+1}}{Y_t} \\log(1+\\pi_{t+1})\n",
    "$$\n",
    "\n",
    "Inflation shows up with two different time displacements, and so we could not express it analytically. Instead, we write a function that returns the residual of the equation, and use the decorator `@solved` to make it into a SolvedBlock."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "@solved(unknowns=['pi'], targets=['nkpc'])\n",
    "def pricing(pi, mc, r, Y, kappap, mup):\n",
    "    nkpc = kappap * (mc - 1/mup) + Y(+1) / Y * np.log(1 + pi(+1)) / (1 + r(+1)) - np.log(1 + pi)\n",
    "    return nkpc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When our routines encounter a solved block in `block_list`, they compute its Jacobian via the the implicit function theorem, as if it was a SHADE model on its own. Given the Jacobian, the rest of the code applies without modification. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Equity price\n",
    "The no arbitrage condition characterizes $(p)$ conditional on $(d, p, r).$\n",
    "\n",
    "$$\n",
    "p_t = \\frac{d_{t+1} + p_{t+1}}{1 + r_{t+1}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "@solved(unknowns=['p'], targets=['equity'])\n",
    "def arbitrage(div, p, r):\n",
    "    equity = div(+1) + p(+1) - p * (1 + r(+1))\n",
    "    return equity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Investment with adjustment costs\n",
    "\n",
    "Sometimes multiple equilibrium conditions can be combined in a self-contained solved block. Investment subject to capital adjustment costs is such a case. In particular, we can use the following four equations to solve for $(K, Q)$ conditional on $(Y, w, r)$.\n",
    " \n",
    " - Production:\n",
    " \n",
    " $$\n",
    " Y_t = Z_t K_{t-1}^\\alpha N_t^{1-\\alpha}\n",
    " $$\n",
    " \n",
    " - Labor demand:\n",
    " \n",
    " $$\n",
    " w_t = (1-\\alpha)\\frac{Y_t}{N_t} mc_t\n",
    " $$\n",
    " \n",
    " - Investment equation:\n",
    "\n",
    "$$\n",
    "Q_t = 1 + \\frac{1}{\\delta \\epsilon_I}\\left(\\frac{K_t-K_{t-1}}{K_{t-1}}\\right)\n",
    "$$\n",
    "\n",
    "- Valuation equation\n",
    "\n",
    "$$\n",
    "(1+r_{t})Q_{t} = \\alpha Z_{t+1} \\left(\\frac{N_{t+1}}{K_t}\\right)^{1-\\alpha} mc_{t+1} - \\left[\\frac{K_{t+1}}{K_t} - (1-\\delta) + \\frac{1}{2\\delta \\epsilon_I}\\left(\\frac{K_{t+1} - K_t}{K_t}\\right)^2\\right] + \\frac{K_{t+1}}{K_t}Q_{t+1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solved blocks that contain multiple simple blocks have to be initialized with the `solved_block.solved` function instead of the decorator `@solved`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "@simple\n",
    "def labor(Y, w, K, Z, alpha):\n",
    "    N = (Y / Z / K(-1) ** alpha) ** (1 / (1 - alpha))\n",
    "    mc = w * N / (1 - alpha) / Y\n",
    "    return N, mc\n",
    "\n",
    "@simple\n",
    "def investment(Q, K, r, N, mc, Z, delta, epsI, alpha):\n",
    "    inv = (K/K(-1) - 1) / (delta * epsI) + 1 - Q\n",
    "    val = alpha * Z(+1) * (N(+1) / K) ** (1-alpha) * mc(+1) - (K(+1)/K -\n",
    "           (1-delta) + (K(+1)/K - 1)**2 / (2*delta*epsI)) + K(+1)/K*Q(+1) - (1 + r(+1))*Q\n",
    "    return inv, val\n",
    "\n",
    "production = solved(block_list=[labor, investment],\n",
    "                    unknowns=['Q', 'K'],\n",
    "                    targets=['inv', 'val'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Making use of these solved blocks, we can write a DAG for this model in just 3 unknowns $(r, w, Y)$ and 3 targets, asset market clearing, fisher equation, wage Phillips curve:\n",
    "\n",
    "![Directed Acyclical Graph](figures/hank2_dag.png) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Determinacy check\n",
    "\n",
    "Let's start by defining the inputs common across our convenience functions. Since computing the Jacobian of this two-asset HA block is somewhat costly, it's a good idea to save it for subsequent use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 300\n",
    "block_list = [two_asset.household_inc, pricing, arbitrage, production, \n",
    "              dividend, taylor, fiscal, finance, wage, union, mkt_clearing]\n",
    "exogenous = ['rstar', 'Z', 'G']\n",
    "unknowns = ['r', 'w', 'Y']\n",
    "targets = ['asset_mkt', 'fisher', 'wnkpc']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are ready to apply the winding number criterion. Recall that a winding number of 0 indicates that the model has a unique solution around the steady state, while a winding number of -1 or less indicates indeterminacy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Winding number: 0\n"
     ]
    }
   ],
   "source": [
    "A = jac.get_H_U(block_list, unknowns, targets, T, ss, asymptotic=True, save=True)\n",
    "wn = det.winding_criterion(A)\n",
    "print(f'Winding number: {wn}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linearized dynamics\n",
    "\n",
    "Computing $G$ is fast using the saved HA Jacobian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = jac.get_G(block_list, exogenous, unknowns, targets, T=T, ss=ss, use_saved=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot some impulse responses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAEYCAYAAAA06gPTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4W/XVwPHvkeVtJ7Zjy9nOcDYBQswOkEGAUmZZCdACb1vaQkop0BI6oNCXUUopu5TSQWkhjAIvFCiFBMIeCYRAFtm2Mx3He9s67x9XTmRHlmVJtuzkfJ5Hj3XX7x5djeN777m/K6qKMcYY09e5Yh2AMcYYEw2W0IwxxuwXLKEZY4zZL1hCM8YYs1+whGaMMWa/YAnNGGPMfsESmjHGmP2CJTRjjDH7BUtoJqpEZJOInBjrOA5E7be9iKwQkemxWHdvbTPIuoJuq970uRaR20Xk6ljH0VNE5GMRmRTKvCEnNBG5VES+EJFaEdkuIn8QkYwuLB+1D0Rv+nD1BuFuDxFJFJE/i8hmEakSkc9E5Gvt5nlLROpFpNr3WBO9yPuGvvp5U9VJqvpWrOPoC9pvq976notIDvAt4I9+47JE5HkRqfF9ly/sYNlOv+9RjDOkmPzmHyEir4hImS+/PCAibt/ku4BbQllvSAlNRK4FfgP8BOgPHAXkAa+LSEIobcSS34YxbbmBIuAEnPf1l8DTIjKi3XzzVDXN9xjXsyH2bfbZM1F2KfCKqtb5jXsQaARygYuAP3SwRxPq9z0aQo2p1UPATmAQcKgvxit8014EZojIoE7XqqpBH0A/oBo4v934NF8A/+MbViDfb/rfgP/1PX8c8AJ1vrZ+CmwCbgBWAmXAX4Ekv+UDtheorQ7i3gRcDywHGnDezMHAv4ASYCNwld/81wNbgCpgDTCrXVsBYwUmAG8B5cAK4IwAcVzni6MCeMpv2YDrDBZngNcZaNsGjamT93s5cI7f8FvAd7qwfLBtFfQ976Ctn/hiqgH+jPMFedW3zd4AMqPwPgT7XATavvOB9b4YVgJnd/LZ+wnwr3bz3A/cE8Y27PB1+pY7MdAwMAx4zvcaS4EHfOO7GltHn9kOt28IcXcUm3/8433vzZzOvq9+7V4GvOQ3vA542m+4CDg0wLo6+r3q8PV1sK3SgUd87+FO4Mehfo+CtLkIuNhvOBUncYxt95m9I5zvezQe4cQErAJO9Rv+LfBHv+HXgUs6XXcIwZ0CNAPuANMeA570Pe8woQX5sn3p+zBnAe+1mz9YgmzTVgdxbwKW+dpPxtkbXQrcCCQAo4ANwMnAON+He7Bv2RHA6M5iBeJ9X5Kf+dqcifMFG9du2Y9xfjSzfG/c9ztaZ7A4O3mtrV/GTmMK0k4uUA+M9xv3Fs4PzS7f654ewnYP+L529p530NaHvriG4PwofApMARJxvtw3Rfg+dLq92feze56vHRdwAU6yHRTkszfIN0+Gb7rb91qmdmUbdvY6A8S5CTgRiAM+B36P82OTBEzzzRNybAT5nnS0fTv7THYSW2v8hwGFwGmdxdEu3lE4CdTle52bgS1+08oAV7Bt19nnp5PvwkKcxJvoe80twMB28/zbF2Ogx78DtFkCHO43PAWoazfPdfgl8q5836PxCCcmnO/i34EUnO/6l/j9owjcB9zd2bpDOeSYDexS1eYA07b5pofrAVUtUtXdwK3A3AjaCuQ+X/t1wOFAjqreoqqNqroB+BMwB+eDlghMFJF4Vd2kqutDiPUonD3VO3xtLsL5gLZ/Hfep6lbfsi/h7FJ3tM5gcYYi1JjaEJF44J/AY6q62m/S9Thf/iE4/22+JCKjO4kh2Pva1ff8flXdoapbgHeAj1T1M1VtAJ7H+fJE8j50eXur6jO+dryq+hSwFjgiwLqKVLVOVbcBb+MkQnD+SdylqkuDvO5IPm/tHYHzQ/wTVa1R1XpVfdf3WroSW2ffk0Dbl07i7jA2n+NwDjldoqr/DjEOfK9tA07ibD2E9RqwRUTG+4bfUVVvJ9vOX0evbx8icpovht+oaoPvNW8BxraL8TRVzejgcVqApjN8r6lVGs4eo78KnL3DDgX5vkdDODEtBiYBlUAxsAR4wW96Fc5rDyqUhLYLyO7gXMAg3/RwFfk934zzwY4m//bzgMEiUt76wPmPMVdV1wFXA78CdorIAhFpH0ugWAcDRe2+FJtxfvz9bfd7XgukBVlnh3GG+JpDjWkPEXHhHBJoBOb5T1PVj1S1yvelfAxnb+HUTmII9r529T3f4fe8LsBwGhG8D4SxvUXkWyKyzG/+g9j3H7uidsOPARf7nl+Ms72DieTz1t4wYHMH/5SGHFsI35NA25dO4u4stu8D76vqm12Iw99iYDpwvO/5WzjJ7ATfcFd09PoCOQP4v9YB33esP20/v+Eoo21iqMY5LeSvH22TXhvBvu9R0qWYfPG8hnPYORXnu5SJU7fRKh1nrzWoUBLaBzjnAb7RLohU4Gs4u9XgvMEpfrMMbNeOBmh7mN/z4cBWv+Fg7QVqKxD/+YqAje3+A0pX1VMBVPUJVZ2G8wOntN2YHcW6FRjme0P8p20JKbjA6wwaZwivs0sxiYiw99zUOara1FnYgHQyT7D3Ndi0cEXyPoSyvfdsXxHJw9mDmwcMUNUMnMMj7bdJ+8/oC8DBInIQcBrOf8fBRPPzVgQMD1KgEnJsIXxPAgkWd2exfd83/fdhxtGa0I7zPV9MaAkt1N+YjhyJcz6w1UycPd82VcIi8qpfBXH7x6sB2l1O2728rwC3iIzxG3cIznnKfYTxfQ9Hl2LCOYQ7DOeoRIOqluKcN/b/Dk7AOTQdVKcJTVUrgJuB+0XkFBGJ91XFPIOza9j639wy4EIRiRORU3A+MP524By68neliAwVkSyc/4qf8psWrL1AbXXmY6BSRK4XkWRfuweJyOEiMk5EZopIIs4x5TqcwxqdxfoRzvmHn/q2y3TgdGBBZ8EEWWeHcQZpzn97dDWmP+B8WE7XtpVTiEiGiJwsIkki4haRi3D+032tk5cX7H0NNi1cYb8PhLa9/bdvKs6PXQmAiFyGs4cWlKrWA88CTwAfq2phJ4tE8/P2Mc7pgTtEJNX3fh7b1dhC/J4EEizuoLHh/Fd/CnC8iNwRRhyLgRlAsqoW4xy2PgUYAHwWJOZwfmPwxRcPjAHO9b2eSThVfNe3n1dVv6Z7K4jbPwKV1L+C32+hqtbg7Nnc4tt+xwJn0vERgA6/79HS1ZhUdRdOwc8PfL8zGcAl+BKY732eilMY0unKQz3R922c/0RbD/v8EV+FmW96AU4GrvIF/iRtizzOxDmxW45zgnATeyu5ynEOe6SE0l77tjqIdxPtCkdwDn08iXPooAyn4OBE4GCcL1YVsBvn+P7gdm0FjBXnuO9inGPEHVW8+Z9c/hXwj2Dr7CjOIO9N+20bNCa/5Vr/u63HOUzQ+rjINz0H+MQXY7kvjtmdfE6Cbaug73ln76Fvu/3Kb/g7wBuRvA+hbO8A2/dW33u2C7jbt97vdLQuv/HTfNv7sgi2YYevM8Br3DOMs0f0As5ewy6c80Fdio3gn9kOt28IcQeMrV38WTg/cL8OFkcHcW8D/uo3vAR4tZPPR6Dfqw5fX4DttNIXawXOOdZLgr3noT5wDscV4yTo1nFZvu1X44v5Qr9prwI/C+X73n7+duvt6vgOYwq0HM75yLdwvn+7cHaYPL5p5wHPhbJ9xLdAjxORTTg/Am/EJIAu6Eux9mYH+nYUkeHAapxKt8og822ih7dTqLGZzonIxTgJ+5xuav82YKeq3tMd7fc2IvIR8G1V/bKzee2iT2N6gO/80TXAgt6WMHpzbH3UIThl/d1CVX/WXW33Rqp6ZKjzWkIzppv5Cqh24FT2nRLjcNrozbH1YQfTeRWr6QYxO+RojDHGRJP1tm+MMWa/ELNDjr5S/Htxur55VFXvaDf9Upz+vFqvsXlAVR8N1mZ2draOGDEi+sEaY8x+bOnSpbtUNSfWcUQqJglNROJwemOejVOC+omIvKiqK9vN+pSqhnwl+4gRI1iyZEkUIzXGmP2fiGyOdQzREKtDjkcA61R1g6o24lxgeWaMYjHGGLMfiFVCG0LbvuqKCdwf3TkislxEnhWRYQGmIyKXi8gSEVlSUlLSHbEaY4zpA2KV0AL1Bdi+3PIlYISqHoxz36vHAjWkqo+oaoGqFuTk9PlDwMYYY8IUq4RWTNvOV4fSrpNaVS1V5xYh4HQGO7WHYjPGGNMHxSqhfQKMEZGRIpKAc++pF/1nkLa32z6Dbrzy3hhjTN8XkypHVW0WkXk4vbbHAX9R1RUicguwRFVfBK4SkTNw7pa9G7g0FrEaY4zpG/arnkIKCgrUyvaNMaZrRGSpqhbEOo5IWU8hwKc7PuXeT++NdRjGGGMiYAkNWFG6gke/eJSSWiv7N8aYvsoSGjAhawIAq3Zb3YkxxvRVltCACQN8Ca3UEpoxxvRVltCAFBI4pnoQK0vbdyVpjDGmr7CEBpQ88CBXPVjMhq0rYh2KMcaYMFlCA9KmHYvLq3hWbqOsvizW4RhjjAmDJTQg+dBD0dRkpqxXO49mjDF9lCU0QOLjST7maKasV1aW2mFHY4zpiyyh+WROn0VWNexY/nGsQzHGGBMGS2g+qcdNAyDh4y9jHIkxxphwWELzifd4qB7pYdSqcioaKmIdjjHGmC6yhOYn7pjDGVcMawo/jXUoxhhjusgSmp/BJ55GnML2t/4b61CMMcZ0kSU0P7mHT6MmWeCDpbEOxRhjTBdZQvMjbjdbJ3kYuHwr6vXGOhxjjDFdYAmtneYjDia9uoXyLz6LdSjGGGO6wBJaOzkzTgKg8L8vxDgSY4wxXWEJrZ3xo49k3SBofPfDWIdijDGmCyyhtZOTksOacamkfFVMc5l1VGyMMX2FJbQAagvG41Koee/9WIdijDEmRJbQAsiacgSVyVCx+M1Yh2KMMSZEltACmJAzic9HCdVvv23l+8YY00dYQgtg0oBJfDpakIoq6r+0zoqNMaYvsIQWQG5KLpvHZaIC1W+/E+twjDHGhMASWgAiwvDhkygelkz122/HOhxjjDEhsITWgQlZE/hwRBP1X3xB8+7dsQ7HGGNMJyyhdWDigIksHaWgSs2778Y6HGOMMZ2whNaBCQMmsHEgNPdPpXqxHXY0xpjezhIaQFMd7FjZZtTQtKGkJfaj+CAPNe++i7a0xCg4Y4wxoYhZQhORU0RkjYisE5H5QeY7V0RURAq6LZh374E/HOMktr3rZWLWRJaMbKGlooK65cu7bfXGGGMiF5OEJiJxwIPA14CJwFwRmRhgvnTgKuCjbg3IMx5QKFnTZvSEARN4zbMTXC6rdjTGmF4uVntoRwDrVHWDqjYCC4AzA8z3a+BOoL5bo/H4cmnJ6jajJ2RNoCKxGSaPo8bOoxljTK8Wq4Q2BCjyGy72jdtDRKYAw1T138EaEpHLRWSJiCwpKSkJL5qsURCXADvbnkebMGACACUHD6N+5Uqaw23fGGNMt4tVQpMA43TPRBEX8Hvg2s4aUtVHVLVAVQtycnLCiyYuHrLHws5VbUbn9csjxZ3C8nw3ANXvWPm+Mcb0VrFKaMXAML/hocBWv+F04CDgLRHZBBwFvNithSE54/dJaC5xMT5rPB+mbcedk2Pn0YwxpheLVUL7BBgjIiNFJAGYA7zYOlFVK1Q1W1VHqOoI4EPgDFVd0m0ReSZARRHUV7YZPXHARNaUfUXKcdOoee89tLm520IwxhgTvpgkNFVtBuYBrwGrgKdVdYWI3CIiZ8Qipr2FIW0rHScOmEh9Sz01BePwVlVRt2xZDIIzxhjTGXesVqyqrwCvtBt3YwfzTu/2gDxOAQg7V8Kww/eMnpDljF+bn8wot5vqxW+TUtB9Rz6NMcaEx3oKaZWRB/Ep+5xHG9F/BElxSXxZv5GUKVOofsduJ2OMMb2RJbRWLhfkjNundN/tcjMuaxyrdq8i7YTjaVi9mqYdO2IUpDHGmI5YQvPnmbjPxdXgHHZcvXs1KccdB2DVjsYY0wtZQvPnmQDVO6CmtM3oiQMmUtNUw/bcBNwDB1Jjd7E2xphexxKav9bCkJK259FaewxZXbaatOOPp+b999HGxp6OzhhjTBCW0PzltFY6tk1oozNGE++KZ2XpStJOOB5vTQ21n1n5vjHG9CaW0Pz1GwyJ/fcpDIl3xTM2cyyrSleRcuRREB9P9duLYxSkMcaYQCyh+RNxDjvuDFAYMmACK3evxJWaQkrBVGqsMMQYY3oVS2jteSY4e2iqbUZPHDCRqsYqiquLSTvueBrWrqNp69YOGjHGGNPTLKG155kA9eVQtb3N6IlZTtdYq0qd69EAqq3a0Rhjeg1LaO35d4HlJz8zH7e4WbV7FQmjRhE/ZIhdj2aMMb2IJbT2Orh7dWJcIvmZ+awqXYWIkHr8cdR8+CFeK983xphewRJae6nZkJqzzx4aOD2GrCxdiaqSdvzxaG0tdUu67442xhhjQhdRQhORY0Uk1ff8YhG5W0TyohNaDHkm7HMtGjiVjmUNZeyo3UHqkUciCQlUL7bDjsYY0xtEuof2B6BWRA4BfgpsBv4ecVSxluMr3fd624xuvZXMytKVuFJSSDn8cOt93xhjeolIE1qzqipwJnCvqt4LpEceVox5JkBTjXMHaz/jssbhEhcrS53DkWknHE/jhg00FhUFasUYY0wPijShVYnIDcDFwMsiEgfERx5WjLUWhrQ77JjsTmZU/1Gs2u2MTzu+tXzfDjsaY0ysRZrQLgAagG+r6nZgCPDbiKOKNc94528HhSGrSp2EljBiBPHDh1vv+8YY0wtEvIeGc6jxHREZCxwKPBl5WD3rhc+2cN7D7+P1+noHSeoP/YYELAyZOGAiJXUllNSWAM5eWs1HH+Gtr+/JkI0xxrQTaUJ7G0gUkSHAQuAy4G+RBtXTahqb+WRTGVsr6vaODFLpCOw97HjCCWh9vR12NMaYGIs0oYmq1gLfAO5X1bOBSZGH1bPyc9IAWLezeu9IzwTY9RW0NLeZd3zWeATZUxiSevRRuAcNonzBgh6L1xhjzL4iTmgicjRwEfCyb1xchG32uDG5TmFm24Q2EVoaoGxjm3lT41PJ65e35zyauN1kXnA+Ne9/QMOGtvMaY4zpOZEmtB8BNwDPq+oKERkFvBl5WD0rKzWBrNQE1pf4JbScIIUhvlvJtMo45xyIj6f8KdtLM8aYWIkooanq26p6hqr+xje8QVWvik5oPSs/J421O/wT2jhAAheGZE1ke812dtfvBsCdk0O/2bMpf+55vLW1PRSxMcYYf9aXo89oTxrrSqrR1vugJaRC5oighSGrS/d2YJx50YV4q6qoePnlfeY3xhjT/Syh+eR70iivbaK0xq/3fM/EoAnN/7Bj8mGHkTh2LGVPPLk3KRpjjOkxltB8xng6qHQsXQfNDW3m7ZfQj6FpQ/dUOgKICJkXXkjDqlXULVvWIzEbY4zZK9Le9kf6eth/TkRebH1EK7ielN9RQtMW2LV2n/knDNjbY0ir/qefhis1lbIn+9y15cYY0+dFuof2ArAJuB/4nd+jzxnUP4nUhLh9Exrsc7NPcHoMKa4upqKhYs84V2oq/c86i6pX/0NzaWl3h2yMMcZPpAmtXlXvU9U3VXVx6yMqkfUwEXEKQ/wT2oAx4HIHLN2fmOV0YLx6d9tkl3nhXLSpifJ/Pdet8RpjjGkr0oR2r4jcJCJHi8hhrY+oRBYD+e0TmjsBBuQHLAwZP8C5Tq39YcfE0aNJOfJIyhcsQFtaujVeY4wxe0Wa0CYD3wXuYO/hxrtCWVBEThGRNSKyTkTmB5j+fRH5QkSWici7IjIxwlg7le9JY3tlPVX1TXtH5owPuIeWlZTFwNSBbSodW2VeeCFNW7fa3ayNMaYHRZrQzgZGqeoJqjrD95jZ2UK++6Y9CHwNmAjMDZCwnlDVyap6KHAncHeEsXYqcJ+OE6FsEzTW7DP/xKyJ++yhAaTPnIHb46HsiSe6K1RjjDHtRJrQPgcywljuCGCdr2eRRmABzl2v91DVSr/BVKDbL+7qsNIRoGTNPvNPGDCBTZWbqG6sbjNe4uPJuOB8at59l8bNm7stXmOMMXtFmtBygdUi8loXy/aHAEV+w8W+cW2IyJUish5nDy1gl1oicrmILBGRJSUlJWG8hL2GZ6WQEOdiXUm7PTTo8N5oAGvK9k12GeeeB243ZQueiigmY4wxoYk0od2Ec9jxNrpWti8Bxu2zB6aqD6rqaOB64BeBGlLVR1S1QFULcnJyQg48EHeci5HZqaz330PLGglxiYErHX0Jzf8C61bxuR7STzyR8ueew1tXt890Y4wx0RVp58SLgdVAuu+xKsSy/WJgmN/wUGBrkPkXAGeFG2dX7FPp6IqDnLEB99Cyk7PJSc4JeB4NnBJ+b0UFla+82l3hGmOM8Ym0p5DzgY+B84DzgY9E5NwQFv0EGOPraSQBmAO0OVQpImP8Br8O7NtdRzcY7UmjcHct9U1+JfeeiQEvrgZfjyG7Aye0lMMPJ3FMPmVPPGH9OxpjTDeL9JDjz4HDVfUSVf0WTrHHLztbSFWbgXnAa8Aq4Gnf/dRuEZEzfLPNE5EVIrIMuAa4JMJYQ5LvScOrsHGXX1WjZwJUboG68n3mnzhgIhsqNlDbtO9tY0SEjLlzqV+xgvovvujOsI0x5oAXaUJzqepOv+HSUNtU1VdUdayqjlbVW33jblTVF33Pf6Sqk1T1UN/lACsijDUkgTsp9hWGBNhLm5A1Aa96O9xL63/GGbhSUih7wvp3NMaY7hRpQvuPr8LxUhG5FHgZeCXysGJnZHYqLmmX0ILcvbpgYAFul5vFRYFPHcalpdHvzDOofOUVmsvKuiNkY4wxRF4U8hPgj8DBwCHAI6p6fTQCi5Wk+DiGZaW0TWj9h0FCGuzcdw+tX0I/jhx4JAsLF3Z4nixz7ly0sZGK56x/R2OM6S5hJzQRiRORN1T1OVW9RlV/rKrPRzO4WMnPaV/p6OqwCyyAWXmzKKwqZG154LqVpLFjSSkooOxJ69/RGGO6S9gJTVVbgFoR6R/FeHqFfE8aG3fV0Nzi3TvSMyFg6T7AjGEzEISFhQs7bDPzogtpKi6m5t13ox2uMcYYonD7GOALEfmziNzX+ohGYLGU70mjscVLUZnfBdGeCVC7C6r37Y0kOzmbKZ4pLNzccUJLnzWLuJxsKw4xxphuEmlCexmnTP9tYKnfo08L2qdjB4cdZw6fyZqyNRRVFQWcLgkJZJ53HtVvv01jcXFU4zXGGBNmQhOR1l2Riar6WPtHFOOLidG+hLZ2Z9XekUFK9wFmDZ8FwKLCRR22m3H++eByUb5gQXQCNcYYs0e4e2iDROQE4AwRmeJ/c8++fIPPVv2S4sntl9h2Dy0tF5IzO9xDG5o+lPFZ44OeR4sfOJD0mTMpf/ZfeBsaoh22McYc0MJNaDcC83H6YLybth0Th3SDz95ujCe9bSfFIpDTcWEIOIcdl+1cxq66XR3Ok3nRhbSUl1P5qvXvaIwx0RRWQlPVZ1X1a8Cdfjf2DPkGn31BvieN9SU1ba8ta6107OB6sxOHn4iiQQ87phx5JAkjR1pxiDHGRFmkF1b/OlqB9DajPWlUNzSzvbJ+70jPBGiohMrANwbIz8hnePrwoIcdRYTMuXOpX76cui++jHbYxhhzwIq0ynG/lZ/jKwzZEdrNPsFJVrPyZvHxto+pbKwMOA9A/7PPQpKTKVtge2nGGBMtltA6EE7pPjjVjs3a3GHfjgBx6en0P/10Kv/9Mi3l+/bgb4wxpusiTmi+LrAGi8jw1kc0Aou17LQEMlLiWVfil9BSspxqxyCFIZOzJ+NJ9gQ9jwbOzT+1oYHy51+IVsjGGHNAi/QGnz8EdgCv41xk/TLw7yjEFXMism+fjuArDOl4D80lLmYMn8G7W96lrrmuw/mSxo8n+bDDKHvySdTr7XA+Y4wxoYl0D+1HwDjffcsm+x4HRyOw3iDfEyihTYSSNRAkCZ2YdyL1LfW8v/X9oO1nXnghTYWF1r+jMcZEQaQJrQioiEYgvVG+J43dNY3srmncO9IzAZrroHxTh8tNzZ1Kv4R+nR527HfSbNyDBrHznnusF35jjIlQpAltA/CWiNwgIte0PqIRWG8QsDAkp7UwpOPzaPGueKYPm86bRW/S5G3qcD5JSMBz7bU0rFxFxfP7xZ13jDEmZiJNaIU4588SgHS/x34hcEIb5/wNch4NnGrHqsYqlmxfEnS+fl8/leQpU9j5+3toqa4OOq8xxpiOuSNZWFVvBhCRdGdQ96tf5MH9k0mOj2ub0JL6Qf/hAe9e7e+YwceQ7E5mYeFCjh58dIfziQi5P7uBTeedT+nDD+O57rpohW+MMQeUSKscDxKRz4AvgRUislREJkUntNhzuYTRntS2ve5D0Jt9tkpyJzFtyDQWFS7Cq8GrGJMnT6b/2Wez+7G/07h5c6RhG2PMASnSQ46PANeoap6q5gHXAn+KPKzeIz8nrW0nxQCe8bDrK2jp+PwYOJ0Vl9SVsLxkeafryfnx1Uh8PDvu/G0k4RpjzAEr0oSWqqpvtg6o6ltAaoRt9ipjctPZWlFPTUPz3pGeieBtgtL1QZc9fujxuF3uTqsdAeI9HgZ873tUL1xIzQcfRBq2McYccCKuchSRX4rICN/jF8DGaATWW4z29em4viRAF1glwQ879kvox5EDj+SNwjfa9trfgaxLLyF+6FB23HY72tzc6fzGGGP2ijSh/Q+QAzwHPO97flmkQfUmrZWObTopzh4L4ur0PBrArLxZFFUVsbZ8bafzuhIT8fz0JzSsXUv5M8+EHbMxxhyIIr19TJmqXqWqh6nqFFX9kaqWRSu43iBvQApul7Tt0zE+GbJGdVq6DzBj2AwEYeHmjm8p4y999mxSjjiCknvvo6Viv71m3Rhjoi6shCYi9/j+viQiL7Z/RDfE2IqPczEyO3XfLrByxoe0h5adnM0Uz5Sg90jz11rG31JZScmDD4YTsjHGHJDCvQ7tcd/fu6IVSG+W70ljzfb2pfsTYc0r0FQP8UlBl585fCZ3LbmLoqoihqUP63R9SeMA16OSAAAgAElEQVTHk3HuuZQ98SSZc+aQOGpUJOEbY8wBIaw9NFVd6nt6qKou9n8Ah0YvvN4h35PG5t21NDT79bfomQDqdcr3OzFr+CyAkKodW+Vc/SNcycnsuOOOLsdrjDEHokiLQi4JMO7SCNvsdfI9abR4lU27aveO7OTu1f6Gpg9lfNZ43tj8RsjrdGdlkX3FFdS8/Q7Vizu+WagxxhhHuOfQ5orIS8DIdufP3gRKQ2zjFBFZIyLrRGR+gOnXiMhKEVkuIgtFJC+cWKOhtXS/zXm0AaPBFR9SYQg4e2mfl3xOSW1JyOvNuuhCEkaMYMcdv0Gbgl/EbYwxB7pw99DeB34HrPb9bX1cC5zS2cIiEgc8CHwNmAjMFZGJ7Wb7DCjw3V/tWeDOMGON2OicNETaJbS4eMgeE9IeGjgJTVHeLHqz85l9JCEBz/U/pXHjRsqeeKKrYRtjzAEl3HNom1X1LVU9ut05tE9VNZQrgo8A1qnqBlVtBBYAZ7Zbx5uq2nqM70NgaDixRkNyQhxDM5Pblu6Dcx6tk4urW+Vn5JPXLy/kasdWadOnkzptGiUPPEjz7t1dWtYYYw4kkXZOfJSIfCIi1SLSKCItIlIZwqJDcG4O2qrYN64j3wZe7SCGy0VkiYgsKSkJ/XBeV+XnBLp79QQoL4SGqsAL+RERZg6fycfbPqayMZRNtHe53PnX462tpeS++7oatjHGHDAiLQp5AJgLrAWSge8A94ewnAQYF7BvKBG5GCgAAvbaq6qPqGqBqhbk5OSEFHQ48j1prC+ppsXrF2ZrYUjJmpDamDV8Fs3azOKirhV5JObnkzl3LuVPP0P9mtDWZYwxB5pIExqqug6IU9UWVf0rMCOExYoB/wuyhgJb288kIicCPwfOUNWGSGONxBhPOo3NXorL/Codc8Y7f0MsDJmcPRlPsqdL5ft7VjXvSuLS09lx+x0h9QtpjDEHmkgTWq2IJADLROROEfkxofW2/wkwRkRG+pafA7TpYUREpgB/xElmOyOMM2KjA929OnMEuJM7vdlnK5e4mDF8Bu9ueZe65rourT8uI4Psq35I7YcfUr2wa+fhjDHmQBBpQvsmEAfMA2pw9rrO6WwhX+HIPOA1YBXwtKquEJFbROQM32y/BdKAZ0RkWay71MoPlNBccZAzLuQ9NIAT806kvqWe97e+3+UYMi+4gMQx+ez4zZ14Gxu7vLwxxuzPwu36CnCqHX1P64Cbu7jsK8Ar7cbd6Pf8xEhii7b+yfHkpCeydp/CkImwPvRDiFNzp9IvoR8LNy/c04NIqMTtxjN/PkXf/g67H3uM7O9+t0vLG2PM/izcC6uf9v39wnfhc5tHdEPsPQJXOo6H6u1QG1pJfbwrnunDpvNW8Vs0ebt+sXTasceSNnMmpX94mOZurOo0xpi+JtxDjj/y/T0NOD3AY780JjeN9Tur2xZleCY5f7d+GnI7s4bPoqqxik+2fxJWHLk//QnepiZ23nNPWMsbY8z+KNwLq7f5nn4DaPZdaL3nEb3wepd8TxpVDc3srPIruBxxLCT2h+Wh35DzmMHHkOxODqvaESBhxAiyvvlNKp57nrovV4TVhjHG7G8iLQrpB/xXRN4RkStFJDcaQfVW+YH6dIxPhoPOhlUvhnSBNUCSO4lpQ6axqHARXvWGFUv2D75PXFYW2375S7z19WG1YYwx+5NI71h9s6pOAq4EBgOLRST0LuX7mNZKx7U72iWuQy6EplpYGXoh5qzhsyipK2F5SXinHOPS0xl06//SsGoV2391s12bZow54EV8YbXPTmA7Tk/7nii12evkpCfSL8m9b5+Ow46ArNHw+ZMht3X80ONxu9xd7tvRX/r06WRfeSUVL7xA+YIFYbdjjDH7g0j7cvyBiLwFLASyge/6esffL4kI+Z4AlY4icMhc2PQOlG0Kqa30hHSOHHQkCwsXRrR3lX3lFaSecDzbb7ud2s8+C7sdY4zp6yLdQ8sDrlbVSap6k6qGfoVxH+UktJp9JxwyBxD4/KmQ25o1fBZFVUV8Vdb5Xa87Ii4XQ+68k/iBA9nyo6tp3rUr7LaMMaYvi/Qc2nwgTUQuAxCRHBEZGZXIeql8Txq7qhsor23XU0fGMBh5nHPYMcQ9rhnDZuASFy+seyGimOL692fo/ffRUlnJlh9fgzaHcgcfY4zZv0R6yPEm4HrgBt+oeOAfkQbVmwXsAqvVIRdC2UYo/DCktrKTszl91Ok8veZpttdsjyiupPHjGXTLzdR+8gk77/pdRG0ZY0xfFOkhx7OBM3D6cURVtwLpkQbVm43xOC8vYEKbcDrEp8Kyf4bc3g8O/QFevDyy/JGIY+t/xhlkXnQRu//2NypfeaXzBYwxZj8SaUJrVKeiQQFEJJSe9vu0IRnJJMW7Aie0xDSYdBaseAEaa/edHqi9tCGcO+Zcnl/7PEWVRZ0v0Inc639K8pQpbP3FL2lYuzbi9owxpq+INKE9LSJ/BDJE5LvAG8CfIg+r93K5hFHZafuW7rc6ZC40VsHql0Nu8/KDL8ftcvPQ5w9FHJ8kJDDknntwpaZQPO+HtFSFdrG3Mcb0dZEWhdwFPAv8CxgH3Kiqodyxuk/L96SxdkcHCS3vWOg/HD5/IuT2clJymDthLi9veJl1Zesiji8+18PQ3/+exi1b2Dr/BtQbXm8kxhjTl0TjjtWvq+pPVPU6VX09GkH1dmM8aWwpr6O2MUA1ocvllPCvfxMqtoTc5v9M+h9S41N5cNmDUYkxpaCA3J/+hOqFCyl9ZL/eaTbGGCD828dUiUhlR49oB9nbtFY6bigJcD0a+K5JU1ge+jVpGUkZfGvit3ij8A1W7IpOh8OZ3/wm/U47jZJ776X63fei0qYxxvRW4fa2n66q/YB7gPnAEGAoTgn//0YvvN4paOk+wIDRMPzoLl2TBvDNid8kIzGD+z+LzlFbEWHQLTeTOGYMW6+9lsbi0PcYjTGmr4n0kOPJqvqQqlapaqWq/gE4JxqB9WZ5A1KJc0nHCQ2c4pBdX8GW0O+TlpaQxrcP+jbvbX2PJduXRCFScKWkMPT++1Cvly1XXWU98xtj9luRJrQWEblIROJExCUiFwEt0QisN0twu8gbkMLanUEqCCedBe6kLhWHAMwZP4ec5Bzu/+z+qPWgn5CXx+A7f0P9ypVsv+XX1jO/MWa/FGlCuxA4H9jhe5znG7ffGxOok2J/Sf1h/GnwxbPQ3NDxfO0XcyfxvYO/x6c7P+W9rdE775U+YwbZV/yAiueeo/ypp6PWrjHG9BaRlu1vUtUzVTVbVXNU9SxV3RSl2Hq1fE8am0traWoJUhJ/6IVQXw5rXu1S298Y8w2GpA3hvk/vi+reVPaVV5J63HFsv/VW6j7/PGrtGmNMbxCt+6EdcPI9aTR7lc2lHVQ6AoyaDumDu3SfNID4uHiuOPQKVu1exRuF0btfqsTFMeS3dxKfm0vxVT+iubQ0am0bY0ysWUILU36O06djhxdYA7ji4ODzYe3rUL2zS+1/feTXGdV/FA989gAt3uidlozLyGDofffSUl5O8VU/wlsTJCEbY0wfYgktTKM9TreVQc+jgXPYUVtgedfOW8W54rjy0CvZULGBlzeG3o1WKJImTmTw7bdR99lnFH77O7RU7veXDhpjDgBRSWgicpSILBKR90TkrGi02dulJLgZkpHccZ+OrXLGwZCpXT7sCHBi3olMyJrAQ8seoqmlKcxIA+t36qkMuef31K1YweZvXWKHH40xfV64PYUMbDfqGpzbyJwC/DrSoPqK/M4qHVsdMhd2fAnblnepfZe4+OGUH7KlegvPr3s+zCg71u+kkxj20EM0btrE5osupmnbtqivwxhjekq4e2gPi8gvRSTJN1yOU65/AXDAHL/K96SxvqQar7eTSsSDzoG4hLD20qYNmcYUzxT++PkfqW+O/kXRacdNY/ifH6V51y42XXQRjZs3R30dxhjTE8Lt+uosYBnwbxH5JnA14AVSgAPikCM4Ca2+ycuW8rrgM6ZkwdhTnPNoXTx0KCJcNeUqdtbt5Kk1ofcN2RUpU6cy/LG/obV1bLr4Yuq/+qpb1mOMMd0p7HNoqvoScDKQATwHrFHV+1S1JFrB9XZjOuvT0d+hF0HtLqfisYsKBhZwzOBjePSLR6luDGFdYUieNIm8fzyOiIvCb36LuuVdOzxqjDGxFu45tDNE5F1gEfAlMAc4W0SeFJHR0QywN+u0k+I2M8+C1Jwud4XV6qopV1HeUM7jqx4Pa/lQJObnk/fEP3Glp1N46WXUfPRxt63LGGOiLdw9tP/F2Ts7B/iNqpar6jXAjcCt0Qqut8tISSA7LSG0hBYXD5PPhzX/gdrdXV7XpOxJzBo+i7+v+Dvl9eVhRBuahKFDyfvnP3APGkTR5ZdTvXhxt63LGGOiKdyEVoGzVzYH2HPFsKquVdU5oTQgIqeIyBoRWSci8wNMP15EPhWRZhE5N8w4u93onLTgnRT7O3QueJuc/h3DMO/QedQ01fCXFX8Ja/lQxefmkvePx0kcPZqiK+dR+WrXuu4yxphYCDehnY1TANJMGJ0Ri0gc8CDwNWAiMFdEJrabrRC4FAjvGF0PGZPrlO6H1OfiwMmQOznsw475mfl8fdTXeXLVk5TUdu+pSndmJsMf+xvJhxzClmuvo/zZ8JKwMcb0lHCrHHep6v2q+rCqhlOmfwSwTlU3qGojsAA4s906NqnqcpzqyV4rPyeNyvpmSqpD7FH/0Ath62ewc3VY67vikCto9jbzpy/+FNbyXRGXns7wR/9E6jHHsO0Xv2T3Y491+zqNMSZcser6aghQ5Ddc7BvXZSJyuYgsEZElJSU9X2A5bmA/ABavCXHdk88DlzvsvbRh/YZx1pizeOarZ9hS3f13oHYlJzP0oQdJP+kkdtx+ByUPPmj3UzPG9EqxSmgSYFxYv5Kq+oiqFqhqQU5OToRhdd2RI7M4eGh/7vrvGmoamjtfIC0H8mfD509BSwjzB/C9g7+HCxcPf/5wWMt3lSshgSF3/47+Z53FrvsfYOedv7WkZozpdWKV0IqBYX7DQ4GtMYolIi6XcNPpk9hR2cBDb60LbaFD50L1dtjwVljrHJg6kAvGX8CL619kQ8WGsNroKnG7GXTbrWRedBG7//pXtt94I9rY2CPrNsaYUMQqoX0CjBGRkSKSgFMt+WKMYonY1LxMzp4yhD+9s5HC0trOFxh7CiRnhn3YEeA7k79DUlwSt390O83e8Pb0ukpcLnJ/8XMGfP97lD/zLJvmXkjDho09sm5jjOlMTBKaqjYD84DXgFXA06q6QkRuEZEzAETkcBEpBs4D/igiK2IRa6iuP2U8bpdw6ysrO5/ZnQgHnQurX4a68K4py0rK4qeH/5QPt33IHR/f0WOHAEUEz9VXM+S+e2kqLmbjN75B2YIFdgjSGBNzMbsfmqq+oqpjVXW0qt7qG3ejqr7oe/6Jqg5V1VRVHaCqk2IVaygG9k/iyhn5vLZiB++t29X5AofOheZ6WBF+L/rnjD2HyyZdxlNrnuLxld3Xg0gg/U46iZEvvkjKYYex/Vc3U3zFlXYLGmNMTNkNPqPo29NGMiwrmZtfWkFzSydXGww+DLLHhdUDv7+rp17N7LzZ3LXkLhYVLoqora6Kz/Uw7NE/kXvDfGree48NZ5xpPYsYY2LGEloUJcXH8fNTJ/LVjmr++VFh8JlFnGvSij6C0vVhr9MlLm6ddisHZR/E/Hfms6K0Z4/MistF1iWXMOKZZ3APGEDR977P9lt+jbeukzsQGGNMlFlCi7KTJ+VybP4A7n79K8pqOqkCPPgCkDh44yaI4BxUsjuZ+2beR2ZiJvMWzmNbdc/fqDNp3FhGPPM0WZdcQtkTT7Dx3POoXxnC+URjjIkSS2hRJiLceNokqhuaufv1Tu4r1m8QzL4ZVr0E79wV0Xqzk7N5cNaD1DfXc+WiK7vtNjPBuBITyb1hPsP+/Cjeqio2XjCH0kcfRb29urMXY8x+whJaNxg3MJ2LjxzOPz/azKptnfQMdvQ8Z09t0a1OT/wRyM/M53fTf8eG8g1c9/Z1PVbO317asccy8v9eIH36dHbe9TsKL72Mpm09v9dojDmwWELrJj+ePZZ+yfHc8tLK4CXtInD6vTDoYHjuu1AS2d2ijxl8DL846he8t+W9Hi3nb8+dmcmQ++5l0K23Uv/ll2w48ywqX3klJrEYYw4MltC6SUZKAtfOHssHG0r5z5fbg88cnwwX/BPiEmDBhVBfEdG6zx17LpcdFJtyfn8iQsY532DkC8+TOHIkW665lq3XX09Ldc8fDjXG7P8soXWjuUcMZ/zAdP735VXUN7UEnzljGJz/dyjbCP/6Lng7mb8TVx+2t5x/YeHCiNqKVMLw4eT98x9kX3klFS/9m41nnkXVwoV2MbYxJqosoXUjd5yLG0+fyJbyOv70dgh9Lo44Fk65A9a+Bm/eFtG6XeLitmm3OeX8b89nxa7YdrQibjc5P5xH3j//gSQkUHzlPDbPvZDaTz6JaVzGmP2HJbRudszobL520EAeems92ypCuDbr8O/AYd9yqh5XvBDRupPcSdw38z6ykrKYtyg25fztpUyZwqiXXmTgLTfTtHUrm7/5LQovv5z6VatiHZoxpo+zhNYDfnbqBFpUuePVEG7qKQKn3gVDj4AXfgDbv4xo3dnJ2Tx04kPUN9dzxcIrYlLO35643WSefz6j//sanuuupW7Z52w8+xtsue4nNBZ2ckG6McZ0wBJaDxiWlcL3jh/F/y3bypJNuztfwJ0IFzwOSf2dIpHaEJYJYnTGaO6efjebKjZx3eLYlfO350pKYsB3vkP+6/9lwOWXU/XGG6w/9etsv+XXNMfgZq3GmL7NEloP+cH00Qzsl8SvXlpBizeEYoj0gXDBP6BqGzxzadg3A2119OCjnXL+re9x+0e396qCjLj+/fFc82NGv/YaGeeeQ9lTT7HupJPZec89tFRVxTo8Y0wfYQmth6QkuLnh1PF8uaWSZ5cWhbbQ0AI47fewcbHTPVaEzhl7Dv9z0P/w9FdP8/eVf4+4vWiLz/Uw6Fe/YvTL/yZ9xnRKH/4j60+cTelf/oq3oSHW4RljejlLaD3ojEMGU5CXyW9fW0NlfVNoC025GI74HnzwAHy+IOIYfnTYj5idN5vfLfkdCzfHtpy/IwkjRjDk7rsZ8a9nSTroIHbeeSfrTz6F8n/9C23uHYdLjTG9jyW0HiQi3HT6JEprGrl/4drQFzz5VhhxHLx4FWz5NKIYWsv5J2dP5trF1/KHZX/oNefU2kueNInhf36U4X/7K26Ph20//wUbzjiT8hdewFtfH+vwjDG9jCW0HjZ5aH/OnzqMv763ifUlIVYcxsXDeX+DtFx46mKo3hlRDEnuJB6e/TCnjjyVhz5/iEtevYTCyt5bXZh61FGMeGoBQ+67F0TYNv8G1p4wnR23307DhhCu7zPGHBCkNxUHRKqgoECXLFkS6zA6VVLVwMy73mLqiEz+dtkRoS+4bTn8+SQYfCh860VwJ0Qcy382/odbPryFZm8z84+Yz9n5ZyMiEbfbXVSV2o8+ouypp6h6YyE0NZFy+OFkzLmA9NmzcSVEvk2MOdCIyFJVLYh1HJGyPbQYyElP5KpZY3hrTQlvru7C3tagg+HMB6DwA/jP9VGJ5ZSRp/DcGc9xcPbB3PT+TVz95tWU1ZdFpe3uICKkHnUUQ3//e8a8uYica66hads2tl57HetOmM7Ou+6ya9mMOUDZHlqMNDZ7OeWetwH4z9XHk+Duwv8Wr98I790Lp90DBZdFJR6venl85ePc++m99E/sz6+P/TXThkyLStvdTb1eat57n/Knn6Jq0ZvQ0kLqMceQccEFpM+cgcTHxzpEY3q1/WUPzRJaDL25eieX/e0TZozL4abTJzEiOzW0Bb0t8MT5sGExfPN5GHlc1GJas3sN89+Zz7rydcwdP5drpl5Dkjspau13t6YdOyn/17OUP/Mszdu2EZeTTcY555B53nnEDxkS6/CM6ZUsofVCfS2hATz6zgbufv0rmluUy6aNYN6MfNKTQtijqCuDP81yeuefeinM/CWkZEUlpoaWBu799F4eX/k4o/qP4o7j7mDCgAlRabunaEsL1W+/TfmCp6h+29kTTj3+ODLO/gZpx03DlRriPw/GHAAsofVCfTGhAeysrOfO19bw7NJistMS+enJ4zh36lBcrk6KM+rK4a3b4eM/QVI/mPkLmHoZuOKiEtcHWz/gF+/+gt0Nu5l36DwunXQpcVFquyc1bd1K+bO+vbaSEiQhgZSjjyJ9xkzSZs4g3uOJdYjGxJQltF6orya0Vp8XlXPzSyv4tLCcg4b046bTJ3H4iBD2unashFd/CpvegdzJcOqdkHdMVGKqaKjg5g9u5vXNrzM1dyq3TbuNwWmDo9J2T9PmZmo//ZTqhYuoWrSIpiKnx5akgw8mfaaT3BLHjOnVVZ7GdAdLaL1QX09o4JSlv/j5Vu54dTXbKuo57eBB3HDqBIZkJHe2IKx8AV77BVQWw+TzYPYt0C/y5KOqvLj+RW776DZc4uLnR/2c00adFnG7saSqNKxdS/WiRVQtepP65csBiB82zElus2aScthhiNsd40iN6X6W0Hqh/SGhtaptbOaPizfw8OL1AHzvhNF8/4RRpCR08gPbWAPv3uNUQbrccPx1cPSVTg/+ESqqKuJn7/yMZSXLOHnEyXz/4O+Tn5kfcbu9QdOOnVS/+SZVixZS+8GHaFMTcf37kzZ9OmkzZ5I27Vg772b2W5bQeqH9KaG12lJexx2vrualz7cyqH8S8782njMOGdz5YbHdG+G1n8OalyFrtHMn7LEnRRxPs7eZP3/xZx5Z/giN3kYOH3g4c8bNYcbwGcS79o/y+JbqGmrefZeqRQupXvw23ooK57zbEUeQUlBAytTDSJo8GVdS36n+NCYYS2i90P6Y0Fp9smk3N7+0gi+3VDI1L5MbT5vIIcMyOl9w3Rvw6nwoXQtjT4GTb4MBoyOOp6y+jOfXPc9Tq59ia81WPMkezht3HueOPZfs5OyI2+8ttLmZ2qWfUr1oIdXvvkfjemePWeLjSTroIFKmHkbyYVNJnnIo7szMGEdrTHgsofVC+3NCA/B6lWeXFnPna2vYVd3AOYcN5coZoxmZnRp8j625ET56GBb/Bloa4ZgfwnHXQkLkh9BavC28s+UdFqxewHtb38PtcjN7+GzmjJ/DFM+U/a7AormsjLrPllH36VJqlyylbsUKaHLunJCQP5qUqQV7klz8kBD2pI3pBSyh9UL7e0JrVVXfxANvruMv726kqUUZ2C+Jo0ZlcdSoARw1agB5A1IC/5BWbYfXb4LlCyB9MEyfD+O/DqnR2aPaVLGJp9Y8xf+t+z+qmqoYlzmOOePncOrIU0mJT4nKOnobb3099V98Qe3ST6n9dCl1n36Gt9rpdNqdm7snuaUcNoXE/HzE+po0vZAltEhXLHIKcC8QBzyqqne0m54I/B2YCpQCF6jqpmBtHigJrdWW8jreXL2TDzeU8uGG3eyqdm6C2WmCK/wIXv0JbPscEBhyGOTPhjGzYfCUiK9jq22q5ZWNr/Dk6if5quwr0uPTOWvMWVww7gLy+uVF1HZvpy0tNKxdS+3SpdQt/ZTapUtp3rHDmeh2kzhqFInjxpE0biyJ48aROHYcbk+O7cmZmLKEFslKReKAr4DZQDHwCTBXVVf6zXMFcLCqfl9E5gBnq+oFwdo90BKaP1VlfUmNL7kFTnBH+hLciAEpiCps+wzWvgHrXofiJYBCygAYPctJbqNnQeqAiGL6bOdnLFi9gNc3v06zNnPs4GOZM34ORw8+msS4yCsveztVpXnrVmqXLaNhzVc0rFlD/Vdf0bxt25554jIyfMlt7N5El5+PK7mTSzWMiRJLaJGsVORo4FeqerJv+AYAVb3db57XfPN8ICJuYDuQo0ECPpATWnuqyoZdNXuS24cbSimpchJcbr9Ejho1gCnDMhiWlcLQzBSGJtaSWvS2k9zWvQG1pTh7b1NhzEkw5kQYNAVc4d2goaS2hGfXPsuza55lZ91OElwJTM6ZTEFuAQUDCzgk5xCS3QfOD3hLRQUNX31F/Z4kt4aGr9aidXXODCIk5OX5Et0YEkeOJH7YcBKGDyOuf//YBm/2O5bQIlmpyLnAKar6Hd/wN4EjVXWe3zxf+uYp9g2v982zq11blwOXAwwfPnzq5s2be+hV9C3BElyrrNQEhmYmMywjkYKEQqY0fMLI8vfpV7ocQSElG/JPdPbeRs0Ia++tydvEB1s/4ONtH7NkxxJW7V6FV724XW4OGnAQBQMLmJo7lSmeKaTGH1jXfanXS1NxMfVr1vjtza2hqbDIuXDex9W/PwnDhhE/bCgJviTXmuzcublImP90mAOXJbRIVipyHnByu4R2hKr+0G+eFb55/BPaEapa2lG7tocWOlWlpLqBLWV1FO951FLk+1tcVkdjsxeATCo53rWckxO/4Fg+p79WAlAXn0F1ynDq++XRkjESskYRnz2a5IFjSM/0EO/u/FxcdWM1n+38jCU7lrBkxxJW7lpJszYTJ3FMyJpAwcACCnILmJI7hX4J/bp1m/RW3tpaGouKaSoqpLGwiMaiQpoKi2gsKqJp61Zobt4zryQkED90qJPwhg93/g4ZjNvjwe3JxZ09AInre/1xmu5lCS2Sldohx17P61V21TS0TXa769iyu4q00i8ZUbOMod5t5MkO8lw7GEwpLtn71lRqCoUMZHvcIErih7A7aSjVqcOpT8vDlZ5LcoKbRLeLpPg4EuNdJLmdvyKNbG1YzabqL1hX+TnrK1fRrE0Iwuj+Yzgst4AJWeMY0T+PvH7DyU7OPqALKrS5maZt22gsLKSpqLhtsissxFtb23YBlwt3djZxubm4c3KI83iI83hw5Xhw5eTgys5BPNsKCTIAAA8LSURBVB5ITQMErypedf4B8iq+YUX3PN87LZx5UAIuE3BZ/Idb52mdv3Va++HWZQMvB85nvX37e9oBvzj3Lqs4z/eNbd/l94yj7WtrjUn94/Yf5zd/2/Ftt1PrtO8dP5oTJ+aG9TmyhBbJSp0E9RUwC9iCUxRyoaqu8JvnSmCyX1HIN1T1/GDtWkLrWfVNLVTWN1FZ10RlVTWNpRvR0g3ElW0ksWozqTWFZNQVk9m0jTi8e5ar1UTKSKdCUynXVCpIpcLvb6XvbzlplJHI7qQqqlN20pxSjCu5CHE17WlLvQlo0wCkaQDSnI20ZONuySHOm4Pbm4E7Lg63y0WcS4gTQcS567WA7zkIvvE4I/dMo928zhw4Pzvs+TFxnvvG7RmmzTDtftAC/Uh5/X74Av0IdvRj3zqP19vuh9jrJb2+mgG1ZWTWVpBVX0FWbSXZ9RVk1VcyoL6CAXWV9Gtql/SA+rh4SpP6szspncqEVCoTUqlITN3zvP1wrTvR2UgRUeLw4kJx4fU9nHHSwTSX7J3H1Wa67jNvmzakbZvity6XKHH4HqKIOH9bx7lEiWtdr6hvudZY2NOGy29+156/Xmec33hpHbdnWPfGL3tja309znT8njvjE4/6NlOmnx3Wlt9fElpMel5V1WYRmQe8hlO2/xdVXSEitwBLVPVF4M/A4yKyDtgNzIlFrKZjSfFxJMXH4UlPAk86jB4EBOjlv6UJygude7ft3kjK7o0k15UyqLYcb1051JUj/9/e2cbKdZx1/Pc/Z99ubNd26oCNE2hTQgVESupGphRaRSqYNIKGRgUFQYkIEqrSClqpFUERVeALlBY+UCjvpQWVYiAEImjUBFrUCpH0xSSpq7SN2zp2EsemjV9i37t3d2cePszsvXv37l7fa3t3b/Y+P2k0M888s+fZOXPm2TPn7EzzCGqepgjN5fUjcDaFc+VmjtZmeKpW40ilwtNl4JnKCZ5pnOC5MtARtEmhYmJHrPPS2GB73MT2sIlNdhmNWGeGGvXYoG51MBFUEk2YCkKOo4mQ4zQsLjo1CTDrcXS2ZDwXBl3dnC8E0uLAJxmFxTRg0h2gyIMtCwOuclmR9dNAunQw6x3MtWRQz7GFHnmkoIGsTsEO6HQozrXgbBvOtuBch+rZNptnZ/nucy9gzYCdjlgz0vO7ZCkFlA0oGlDWu7FRVo2iGimqkbJqlNVAUTPKskNZixSVQFkJFIS1dr/1gQ1Jr4hAxfKwIB9Sfj6dKy75t3vR4X+sdtYX7SY0k5NL8cmedI7bs9CZXwwhxaHT5ERocsSaHLE2RxU4UkSOFuJopWBuwMsSMuMlMbItRraGyPYY2RoC22NkW4hsjYHtIbI1Rl4SIjNmNMyoW2QmGlW6bmuCqACViwNbUfYNemWfvExlS/Ldcg3RLzCJ2BKhaYSm0Zk1wlwgzEXCbIfOXCDMBsJshzDbTnGzA+E8Y4ygqFcpLqtRztQoZlKsRi3J61WKeo2iUU+yRkoX9Rpq1ClmGimf08p5VarJbgqiRMQIQCxEQERIeUQAAjHllWSp3HriXB9blFvKL8gsLuhEs4V6C2mLBGKKYyBaJA7IL6RzCLY8b9iSOrd+763s3bX3wrqQ36E5zgioNqC6E7bsXHPVEtiVww/1lZkZ3zr7HCfOPcup+VOcap7kVPMUp1qnODl/itPzpznVOsPx+dN8tXWGU60zNMP88oP0IUSjrDFT1qkXNRplnUbZjZNspqxRL+uUKqkUJWVRUqikUlQoVVKqQlmUlEWFsqhQUUmRZZWiQqGSIjsXLfxKF+o6pq4t0oJNvfk0LWlL4kjM052WpzvjMp2ufEkgLpcNC0RiDDDfppydp5ydpzI7T2W2RXW2RWWuRXW2TXWuRW2uQ3WuQ22uTX1ultrJF6jOB6qtSLUVqbUilfM5xj6CoF2BVqUvLlPcrmhRXi7qdMqeUGghHYqlZf35TiFCAaFMTjEUKcSiL63lcuvOa2cKFRQqKFUupAfle2U3XnnjmtpnGnGH5mwIJHHFll1csWXXqus0O83k/HI4M3+G+TDPXGeO+TBPs9OkGZop7kl3y86FOZ5vn6HZTPLur+xO7BAsEGJIcZa/mChUUFAsG1wlLQywIqeLPOhWCsqtJcW23oF4hlKbF3U0PK6YqLeMWitSb0GtndK1VqQ6H6m2QnJ884FKO1J24kJctgNlO1JtB7a0A0U7ULYCRbtDcbZD0Q6o3aGYb6NOQCGiTmAN84gXT1mmv1wURUpLS2KKIqcLKLSgq0JQiCt+FXj5+Mxdj7hDc5whNCoNdlZ2snPT2u8W10rvtFKIgY51Fh1eDAsvonRZfAmlL+6Tdx/ldR1MegEmxV0H1JUVFAt3dV39Qc6qqz/tmBm021i7jXU6Ke5Pt9pYp72Yb7chRqzTgRCwToCYYguLsmXpTsBCSHVjyK9hBixa0rPYI4sQIljEQsx1ov/hHndojrMu6DoNIM2dOhNHEtRqvqD0iwhfUsBxHMeZCtyhOY7jOFOBOzTHcRxnKnCH5jiO40wF7tAcx3GcqcAdmuM4jjMVuENzHMdxpgJ3aI7jOM5UMFWLE0v6P+BCt6zeAXzrvFrjx+1aG27X2lmvtrlda+Ni7PoeM3vRr9c/VQ7tYpD0hfW42rTbtTbcrrWzXm1zu9bGerVrnPiUo+M4jjMVuENzHMdxpgJ3aIv8+aQNGILbtTbcrrWzXm1zu9bGerVrbPgzNMdxHGcq8Ds0x3EcZypwh+Y4juNMBRvOoUm6SdJXJR2SdNeA8rqk/bn8EUkvG4NNV0n6tKQnJH1Z0q8N0LlR0mlJj+bw3lHblY97WNKX8jG/MKBckv4wt9fjkvaMwaZX9rTDo5LOSHpnn87Y2kvShyWdkHSwR3a5pIckPZnj7UPq3p51npR0+4hter+kr+TzdJ+kbUPqrnjOR2TbPZKe6TlfNw+pu+L1OwK79vfYdFjSo0PqjqTNho0Nk+5f6xYz2zCBtBfw14GrgRrwGPADfTp3An+a07cB+8dg1y5gT05vAb42wK4bgX+bQJsdBnasUH4z8AAg4DXAIxM4p8+R/hg6kfYCXg/sAQ72yH4PuCun7wLeN6De5cA3crw9p7eP0KZ9QCWn3zfIptWc8xHZdg/w7lWc6xWv30ttV1/57wPvHWebDRsbJt2/1mvYaHdoe4FDZvYNM2sBfw/c0qdzC/DRnP4n4A2SNEqjzOyYmR3I6ReAJ4DdozzmJeQW4G8s8TCwTdKuMR7/DcDXzexCV4i5aMzsM8DzfeLefvRR4KcHVP0J4CEze97MTgIPATeNyiYze9DMOjn7MHDlpTjWWhnSXqthNdfvSOzKY8DPAh+/VMdbpU3DxoaJ9q/1ykZzaLuBoz35p1nuOBZ08sV/GnjpWKwD8hTnq4BHBhT/sKTHJD0g6QfHZJIBD0r6oqRfGVC+mjYdJbcxfJCZRHt1+U4zOwZpUAK+Y4DOJNvuDtKd9SDOd85HxTvydOiHh0yhTbK9XgccN7Mnh5SPvM36xob13r8mwkZzaIPutPr/t7AanZEgaTNwL/BOMzvTV3yANK12HfBB4F/GYRPwI2a2B3gj8HZJr+8rn2R71YA3Af84oHhS7bUWJtJ2ku4GOsDHhqic75yPgj8BXgFcDxwjTe/1M7G+BvwcK9+djbTNzjM2DK02QDbV/9PaaA7taeCqnvyVwLPDdCRVgK1c2PTImpBUJXXYj5nZP/eXm9kZMzub058AqpJ2jNouM3s2xyeA+0jTPr2spk1HxRuBA2Z2vL9gUu3Vw/Hu1GuOTwzQGXvb5RcDfhL4ecsPWvpZxTm/5JjZcTMLZhaBvxhyzIn0tTwO3ArsH6YzyjYbMjasy/41aTaaQ/s8cI2kl+df97cB9/fp3A903wZ6C/CpYRf+pSLPz/8V8ISZ/cEQnZ3dZ3mS9pLO3bdHbNcmSVu6adJLBQf71O4HflGJ1wCnu1MhY2Dor+ZJtFcfvf3oduBfB+h8EtgnaXueYtuXZSNB0k3ArwNvMrPZITqrOeejsK33ueubhxxzNdfvKPgx4Ctm9vSgwlG22Qpjw7rrX+uCSb+VMu5Aeivva6S3pe7Ost8mXeQADdIU1iHgc8DVY7DpR0lTAY8Dj+ZwM/A24G1Z5x3Al0lvdj0MvHYMdl2dj/dYPna3vXrtEvDHuT2/BNwwpvN4GclBbe2RTaS9SE71GNAm/Sr+ZdJz1/8Enszx5Vn3BuAve+rekfvaIeCXRmzTIdIzlW4f677N+13AJ1Y652Nor7/N/edx0mC9q9+2nF92/Y7Sriz/SLdf9eiOpc1WGBsm2r/Wa/ClrxzHcZypYKNNOTqO4zhTijs0x3EcZypwh+Y4juNMBe7QHMdxnKnAHZrjOI4zFbhDc5wxImmbpDsnbYfjTCPu0BxnTEgqgW2kHR3WUk+S/Fp1nPPgF4njDEHS3Xnvrf+Q9HFJ75b0X5JuyOU7JB3O6ZdJ+qykAzm8NstvzPtZ/R3pj8O/C7wi75v1/qzzHkmfzwvz/lbP5z0h6UOkdSmvkvQRSQeV9t161/hbxHHWN5VJG+A46xFJryYtrfQq0nVyAPjiClVOAD9uZk1J15BWnbghl+0FrjWzb+YV0681s+vzcfYB12QdAffnhW2PAK8kre5wZ7Znt5ldm+sN3JzTcTYy7tAcZzCvA+6zvOahpPOtGVgF/kjS9UAAvq+n7HNm9s0h9fbl8L85v5nk4I4AT1naYw7S5oxXS/og8O/Ag2v8Po4z9bhDc5zhDFoXrsPiVH2jR/4u4DhwXS5v9pSdW+EYAn7HzP5siTDdyS3UM7OTkq4jbdr4dtJmk3es5ks4zkbBn6E5zmA+A7xZ0kxeSf2nsvww8OqcfkuP/lbgmKXtT94KlEM+9wVgS0/+k8Adeb8rJO2WtGyzxrz1TWFm9wK/Cey5oG/lOFOM36E5zgDM7ICk/aTVzZ8CPpuLPgD8g6S3Ap/qqfIh4F5JPwN8miF3ZWb2bUn/Lekg8ICZvUfS9wP/k3e7OQv8AmnaspfdwF/3vO34Gxf9JR1nyvDV9h1nFUi6BzhrZh+YtC2O4wzGpxwdx3GcqcDv0BzHcZypwO/QHMdxnKnAHZrjOI4zFbhDcxzHcaYCd2iO4zjOVOAOzXEcx5kK/h9G+R+qW7Or7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rhos = np.array([0.2, 0.4, 0.6, 0.8])\n",
    "drstar = -0.0025 * rhos ** (np.arange(T)[:, np.newaxis])\n",
    "dY = 100 * G['Y']['rstar'] @ drstar\n",
    "\n",
    "plt.plot(dY[:21])\n",
    "plt.title(r'Output response to 25 bp monetary policy shocks with $\\rho=(0.2 ... 0.8)$')\n",
    "plt.xlabel('quarters')\n",
    "plt.ylabel('% deviation from ss')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear impulse responses\n",
    "\n",
    "Let's compute the nonlinear impulse response for the $\\rho=0.6$ shock above. (Don't forget to use the saved Jacobian!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "On iteration 0\n",
      "   max error for asset_mkt is 4.22E-06\n",
      "   max error for fisher is 2.50E-03\n",
      "   max error for wnkpc is 6.15E-08\n",
      "On iteration 1\n",
      "   max error for asset_mkt is 2.79E-04\n",
      "   max error for fisher is 1.49E-06\n",
      "   max error for wnkpc is 2.61E-05\n",
      "On iteration 2\n",
      "   max error for asset_mkt is 8.78E-06\n",
      "   max error for fisher is 9.99E-08\n",
      "   max error for wnkpc is 7.67E-07\n",
      "On iteration 3\n",
      "   max error for asset_mkt is 5.15E-07\n",
      "   max error for fisher is 2.62E-09\n",
      "   max error for wnkpc is 2.08E-08\n",
      "On iteration 4\n",
      "   max error for asset_mkt is 3.12E-08\n",
      "   max error for fisher is 1.39E-10\n",
      "   max error for wnkpc is 1.03E-09\n",
      "On iteration 5\n",
      "   max error for asset_mkt is 1.92E-09\n",
      "   max error for fisher is 7.90E-12\n",
      "   max error for wnkpc is 5.64E-11\n"
     ]
    }
   ],
   "source": [
    "td_nonlin = nonlinear.td_solve(ss, block_list, unknowns, targets,\n",
    "                               rstar=ss['r']+drstar[:,2], use_saved=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see rapid convergence and mild nonlinearities in the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecFPX5wPHPs3u9Uu6oR5Ui7eiIBcHY0Cj2rrEkGn+KJpqY6E/FEvOLJTGWGFtiMJrYMEZUjB1LVAQU6QrCASdF6lEOruw+vz++c8dyXtk7bm/ubp/367Wv3Zn5zswzs7Pz7Hxn5juiqhhjjDEAAb8DMMYY03xYUjDGGFPJkoIxxphKlhSMMcZUsqRgjDGmkiUFY4wxlSwptCIi8r8i8he/4zCmuRCRqSJyh/d5nIh81YTzVhHp08jTrFyeWGnVSUFEzhWROSKyU0TWicjrInKY33E1BhGZICKFkf1U9f9U9Sd+xdQaiMitIvL0fox/hIi8JyJFIlJQZViCiDwrItu8bTEzYtiNInLNfoTeZERkpoi0uO1MVT9U1f5+x9HctdqkICLXAvcB/wd0BLoDfwZO8jOu1kBEEvyOoRnbBTwBXFfNsFMBBXKA7cBPAUSkF3Ai8GATxegr236aOVVtdS8gG9gJnFFLmWRc0ljrve4Dkr1hE4BC4BfAd8A64OKIcY8HFgM7gG+BX3r9LwI+qjIfBfp4n6fiEtPrXnz/BTp5894KLAWGR4xbANzgzWsr8DcgBUgHdgNhbzo7gS7ArcDTEeNPAhYB24CZwIAq0/4lMB8oAp4DUmpYVxd5sf4R2ALc4fW/BFjixfYG0MPrL17Z77xpzwcGR6yDR4C3vPX3fsV43vBDgNneeLOBQyKGzQR+48WyA3gTyPGGpQBPA5u95Z0NdIzYHv7qfY/fAncAwWqWcyJQCpR56/RLr38XYLq37MuBS6PYBo8CCqr0+zXwU+/z5cCfvc+vAIdFMc2ZXuwfe/G9ArQH/oFLMrOBnvu7Lr3hY735bAO+BCZ4/X8LhIA9Xgx/8vrfD6zx4pgLjIuY1q3ANO/72Q7cBBQD7SPKjAQ2AonVLHfF+M95sX4ODI0YPsBbnm247X1SxLCp7N1eJwCFEcO6Af/y5rsZ+BNuv7AFGBJRrgPu95ZbTWx9cNtwEbAJeK7Kb/9yYBnuN/IQIN6wgLceVuF+J38HsiPGPSxi/a8BLqpmeTKB94AHKqbbKPvPxppQc3rhftzlQEItZW4HPvW+8FzvC/hNxMZT7pVJxCWBYqCtN3xdxUYPtAVGeJ8vou6ksMn7AaQA7wIrgR8BQdwP/r2IcQuAhd7G2w73A652A4/48Tztfe6H+9d6tLcMv8Lt0JIipv0ZbofXDrdzv7yGdXWRtz6uAhKAVOBkb3oDvH43AR975Y/F7Rja4BLEAKBzxDrYARyO+wHeX7HOvDi2Ahd40zzH627vDZ8JfOMtW6rXfac37Ke4nWSaty5HAlnesH8Dj+KSaQdvuX9aw7JWrsOIfu/jknkKMAy3Ezmyjm2wuqTwQ9yOLcl7vxI4BfhblNv1TG+dH4BLdIuBr715JeB2LH9rhHXZFbeTPB638zra686NGPcnVWI7H5egEnB/ptbj/cnw1mmZt80EvPnNAP4nYvw/Ag/W8p2UAafjtuVf4n43id5rOfC/3nr9AW776h+xvX3vN+NtI1968033vtvDvGF/Bu6KmP/PgFdqiO0Z4EZvuSqnEfHbfxX3O+jubTcTvWGXeHH3BjJwyekpb1h3bxnO8ZavPTAscnm8fp9VLFuj7j8be4LN4QWcB6yvo8w3wPER3cfi/Yi9jWc3EUkFl83Hep9X43ZCWVWmeRF1J4XHI4ZdBSyJ6B4CbIvoLiBiR437kX5TdQOv8uOpSAo3A89HDAvg/iVPiJj2+RHD7wYeqWFdXQSsrtLvdeDHVaZfDPTA/TC/xv3bDFQZbyrwbER3Bu6fZzfcDuyzKuU/Ye+/pJnATRHDrgD+432+BJfY86uM3xEoAVIj+p1DRPKtaR163d28+DIj+v0OmFrH9lVdUhDgTtyR02O4H/Y8XKL6LfABboeUVMM0ZwI3RnT/AXg9ovtEYJ73eX/W5a/xdlARw98ALowY9yd1LP9WvH/z3jr9oMrws4D/ep+DuCQyppbv5NMq29o6YJz3Wh+5neF21LdGbG/VJYWDcTvp7/1xBA7C/TsPeN1zgDNriO3v3neZV80wZd8k8Txwvff5HeCKiGH9cYkvAVc78FIN85uKq55cCFxX23fQ0FdrPaewGcipo+6yC+7QrcIqr1/lNFS1PKK7GLcDAzgNt4NeJSLvi8jB9YhtQ8Tn3dV0Z+xbnDW1xFibfZZPVcPetLpGlFkf8Tly+aqzpkp3D+B+76TpNtwhtwBdVfVd3KH4Q8AGEXlMRLKqm5aq7vTG7VI1Zs+qKGN+CrfjelZE1orI3SKS6MWZCKyLiPVR3I44Gl2ALaq6o5aYoqLO9aqar6qXAdfjqtJGea/xuH+7l9QymWi3n/1Zlz2AMyrWl7fODgM61xSUiPxCRJZ4J9i34Y5kciKKVN1+XgYGikhv3JFIkap+VtP02XebCeOqdyu2mTVev5qWszrdgFVVfuMV05+FO8oeLyIH4qqIptcwnV/htvvPRGSRiFT97mpax9XtfxJwf2K64f601uSHuKOtR2op02CtNSl8gqvzPLmWMmtxG3+F7l6/OqnqbFU9Cbdj+TfuHwC4DSmtopyIdKpHzDXpVkOMWsd4+yyfiIg3rW8bGEfV+a3BVcG0iXilqurHAKr6gKqOBAbhqigiT7xWLpOIZOCqOirO7UR+J+CWuc6YVbVMVW9T1YG4uvQTcNVya3BHCjkRcWap6qAol3Mt0C7ySqFoY6qNiAz24nwMd4Q4V91fwdlA/v5M29PgdYlbZ09V+W7TVfVOb/g+60hExuGOLs7EVbG2wdWxS0SxfcZR1T243815uKOap+qIKXKbCQB57N1munn96rOca4DutfxxfBJXJXYBMM2L93tUdb2qXqqqXXC1B3+O8jLU6vY/5bgkvwZXRViTx4H/ADNEJD2KedVLq0wKqloETAEeEpGTRSRNRBJF5DgRudsr9gxwk4jkikiOV77OSxFFJElEzhORbFUtw504C3mDvwQGicgwEUnBHfburytFJE9E2uHqTZ/z+m8A2otIdg3jPQ/8UESO9P4x/wK3c/y4EWIC9y/lBhEZBCAi2SJyhvd5tIgc5M13Fy5BhyLGPV5EDhORJNzJzlmqugZXz9zPu5Q4QUTOAgbi6mVr5V0KOkREgrjvpAwIqeo63EnUP4hIlogEROQAERlfw6Q2AD0rdjJeXB8DvxORFBHJB36MO7lbXRwB77tPdJ2S4i1nZBnBHUX9zPuHuxKoWB/jgRV1LW8UGrwucb+DE0XkWBEJesswQUTyvOEbcHXhFTJxO7SNQIKITAGyqNvfcVWTk6j7tzdSRE71duI/x23LnwIV/+p/5f3GJ+Cq0Z6tY3qf4aqg7hSRdG8ZD40Y/hTufM/5XpzVEpEzItbLVlzyC9VUPsIzwDUi0sv7Y/R/uJPU5bht6ygROdP77tqLyLAq408GvgJeFZHUKOYXtVaZFABU9V7gWtwJ0I247DsZ988e3MmaObj63QW4KxqivSnkAqBARLbjri4435vn17iT02/jrjj4qBEW5Z+4ndoK73WHN6+luA1rhXeIv0+1kqp+5cX1IO7k9onAiapa2ggxoaovAXfhqmu24+o4j/MGZ+H+zWzFHRZvBn5fZZluwVUbjcT9W0RVN+P+4f/CG+dXwAmquimKkDrhrlDZjjtp/j57dzQ/wlXLVFzFNY2aq0Je8N43i8jn3udzgJ64f3cvAbeo6ls1jH84rhpnBu7f327c9xfpYmChqs7xuv/lTXsj7jzDo7Uvat32Z116ifAk3J+Qit/OdezdX9wPnC4iW0XkAVy13eu480ircH8CqlYXVTef/+KuoPtcVQvqKP4y7jxExcnzU72jw1JcUjkOt53/GfiR9/uobd4h3G+iD+4cYaE3/Yrhhbh9ggIf1jKp0cAsEdmJq2L6maqurGNZwJ0XeAp3Hmklbp1d5c17Na56+he438g8YGiV+BW4DLeeX/b+iDSKisujTDMk7uann6jq237H0lhEZCruZN9Nfsdi/Cci7wL/VNUa78QXkVtxF2uc32SBufk+AayNt23VbiIxxvhCREYDI2iGN5SKSE/czYbD/Y2k6bXa6iNjTPMlIk/iqll/XuXKLt+JyG9w1aH3RFkV1KpY9ZExxphKdqRgjDGmUos7p5CTk6M9e/b0OwxjjGlR5s6du0lVc+sq1+KSQs+ePZkzZ07dBY0xxlQSkap3uFfLqo+MMcZUsqRgjDGmkiUFY4wxlVrcOQVjTOtUVlZGYWEhe/ZU2/aciVJKSgp5eXkkJiY2aHxLCsaYZqGwsJDMzEx69uyJazPQ1JeqsnnzZgoLC+nVq1eDphHT6iMRmSgiX4nIchG5vprhF4nIRhGZ571a3MPAjTGNY8+ePbRv394Swn4QEdq3b79fR1sxO1LwmjB+CPcAjUJgtohMV9XFVYo+p6qTYxWHMablsISw//Z3HcbySGEMsFxVV3jN2z5LM2z4yhhjzF6xTApd2bdN9UKqf0TeaSIyX0SmiUi3aoYjIpeJyBwRmbNx48aGRfPGjfDoeHjtFw0b3xjT6mVkuKdlrl27ltNPP93naPwRy6RQ3TFM1db3XgF6qmo+rsXEJ6ubkKo+pqqjVHVUbm6dd2lXb92XsG4eFM5u2PjGmLjRpUsXpk2bFtN5lJd/7/HQzUIsk0Ih+z5fuOKZqpVUdbOqlnidj+OewhUbnb0HF323BEJlMZuNMablKygoYPDgwQBMnTqVU089lYkTJ9K3b19+9atfVZZ78803OfjggxkxYgRnnHEGO3fuBOD2229n9OjRDB48mMsuu4yK1qgnTJjA//7v/zJ+/Hjuv//+pl+wKMTyktTZQF8R6YV7iPbZwLmRBUSks/cMXXCP1FsSs2g6DXHvoVLYuHRvtzGm2bntlUUsXru90ac7sEsWt5w4qN7jzZs3jy+++ILk5GT69+/PVVddRWpqKnfccQdvv/026enp3HXXXdx7771MmTKFyZMnM2XKFAAuuOACXn31VU488UQAtm3bxvvvv9+oy9WYYpYUVLVcRCbjnt8aBJ5Q1UUicjswR1WnA1eLyCTcQ7+34B7iHRud8vd+XjffkoIxzdjitduZtXKL32FUOvLII8nOzgZg4MCBrFq1im3btrF48WIOPfRQAEpLSzn44IMBeO+997j77rspLi5my5YtDBo0qDIpnHXWWdXPpJmI6c1rqjoD9wDzyH5TIj7fANwQyxgAykJh5hfnMiyQRDBcCusXxHqWxpj9MLBLVrOabnJycuXnYDBIeXk5qsrRRx/NM888s0/ZPXv2cMUVVzBnzhy6devGrbfeus99A+np6Q0LvonExR3NxSUhTnv0M15OymNoYAWsn+93SMaYWjSkiqepjR07liuvvJLly5fTp08fiouLKSwspEOHDgDk5OSwc+dOpk2b1qKuZIqLpJCdlkiP9mksKurhJYUFEA5DwNoDNMY0TG5uLlOnTuWcc86hpMRdL3PHHXfQr18/Lr30UoYMGULPnj0ZPXq0z5HWT4t7RvOoUaO0IQ/ZmfzPz2mz6O/ckfg31+PqL6Bd70aOzhjTUEuWLGHAgAF+h9EqVLcuRWSuqo6qa9y4+aucn5fNR+HB3FJ2IVvOfBkyu/gdkjHGNDtxUX0EMKRrGwq0MwWhzoxjAEclpvgdkjHGNDtxc6QwuOveqw7mf1vkYyTGGNN8xU1SyExJpHeuuxRsQeE2n6MxxpjmKW6SAkB+12xGyVIuWn0Deu8g2NnAxvWMMaaViqukMCSvDalSynidg2wvhPVf+h2SMcY0K3GVFPLzslkU7rm3h93ZbIyJgZkzZ3LCCScAMH36dO68806fI4pe3Fx9BDCwcxZbJYt12o7OssW1gWSMMTE0adIkJk2aFNN5hEIhgsFgo0wrro4U0pMT6JObweJwD9fDmrswxkQoKChgwIABXHrppQwaNIhjjjmG3bt3M2/ePMaOHUt+fj6nnHIKW7duBVxT2L/+9a8ZM2YM/fr148MPP/zeNKdOncrkye6JwxdddBFXX301hxxyCL17997nmQ333HMPo0ePJj8/n1tuuaWy/8knn8zIkSMZNGgQjz32WGX/jIwMpkyZwkEHHcQnn3zSaOsgro4UAIbkZbNoSw+O5At08zdIyU5IzvA7LGNMVV/8A+b9s/YynYbAcRFVM+vmw39qaGNz2Lkw/Lw6Z7ts2TKeeeYZHn/8cc4880xefPFF7r77bh588EHGjx/PlClTuO2227jvvvsA97Cczz77jBkzZnDbbbfx9ttv1zr9devW8dFHH7F06VImTZrE6aefzptvvsmyZcv47LPPUFUmTZrEBx98wOGHH84TTzxBu3bt2L17N6NHj+a0006jffv27Nq1i8GDB3P77bfXuUz1EXdJIb9rNp/M6wmAoLBhEXQ/yN+gjDHft201rPqofuPsKap5nJ6HRTWJXr16MWzYMABGjhzJN998w7Zt2xg/fjwAF154IWeccUZl+VNPPbWybEFBQZ3TP/nkkwkEAgwcOJANGzYA7mE9b775JsOHDwdg586dLFu2jMMPP5wHHniAl156CYA1a9awbNky2rdvTzAY5LTTTotqmeoj7pLCkLw2/EV77u2xfr4lBWOaozbdoUcdO/Kqz0VJya55nDbdo5pt1Wayt22r/b6mivIVTWrXZ/oVbc+pKjfccAM//elP9yk7c+ZM3n77bT755BPS0tKYMGFCZTPcKSkpjXYeIVLcJYWBnbNYJx3YrmlkSbF7drMxpvkZfl5U1T376JwPF7/WqGFkZ2fTtm1bPvzwQ8aNG8dTTz1VedTQWI499lhuvvlmzjvvPDIyMvj2229JTEykqKiItm3bkpaWxtKlS/n0008bdb7VibukkJoUpG+HTO7eeBadOnZk8riz/Q7JGNPMPfnkk1x++eUUFxfTu3dv/va3vzXq9I855hiWLFlS+eS2jIwMnn76aSZOnMgjjzxCfn4+/fv3Z+zYsY063+rETdPZkX417Uuen1NIZkoC8285BhFppOiMMQ1lTWc3Hms6u56G5LUBYMeeclZtLvY5GmOMaT7iMinkd82u/GwtphpjzF5xmRQO7JxJYlC4KeEpRrx5BrzTuNf5GmMapqVVZzdH+7sO4zIpJCcE6d8pk4MCS8jbtRBWz/I7JGPiXkpKCps3b7bEsB9Ulc2bN5OS0vCHiMXd1UcVhnRtw6INPRkSKEDXz0dUwU44G+ObvLw8CgsL2bjRmrTfHykpKeTl5TV4/LhNCvl52Sya2xMAKdkO21ZB256+xmRMPEtMTKRXr15+hxH34rL6CGBI1yrNaFuLqcYYE79JoV/HTFYEexJWr8rInq1gjDHxmxSSEgL06NyBldrJ9bBmtI0xJn6TArj7FRare7aCWvWRMcbEd1IYEvF4TtmxFnZt8jcgY4zxWUyTgohMFJGvRGS5iFxfS7nTRURFpM52ORpTfl4274eH8tuyc3l/7F8gyR62Y4yJbzG7JFVEgsBDwNFAITBbRKar6uIq5TKBq4Emv4OsT24GKxN6saSsByUlPRif2PAbPowxpjWI5ZHCGGC5qq5Q1VLgWeCkasr9Brgb2BPDWKqVEAwwqItrB2l+obWBZIwxsUwKXYE1Ed2FXr9KIjIc6Kaqr9Y2IRG5TETmiMicxr7bcYjXON7iddspC4UbddrGGNPSxDIpVNdmRGWjJiISAP4I/KKuCanqY6o6SlVH5ebmNmKI7rzCQbKEP8q9hB8YAaW7GnX6xhjTksQyKRQC3SK684C1Ed2ZwGBgpogUAGOB6X6cbG4rO/hh8DOSi1bChkVNOXtjjGlWYpkUZgN9RaSXiCQBZwPTKwaqapGq5qhqT1XtCXwKTFLV/XusWj31ysmgICGivRV7ZrMxJo7FLCmoajkwGXgDWAI8r6qLROR2EZkUq/nWVzAgZHfpy3ZNdT2suQtjTByLaSupqjoDmFGl35Qayk6IZSy1GZLXliVre3CQLCW8bn5839FnjIlrtv9j3zub+W4xhMp9jccYY/xiSQHIz2tT2QZSIFQCm772OSJjjPGHJQWgR7s0ViYesLeHtZhqjIlTlhSAQEBI6zKIEvVOsdjJZmNMnIrbx3FWNbBbe/6w6gy2SjZ3DLuIZL8DMsYYH1hS8OR3bcOVoRMBOKckhxE+x2OMMX6w6iNPfl525ecF1jieMSZO1ZkURORQEUn3Pp8vIveKSI/Yh9a08tqm0iYtEbAWU40x8SuaI4WHgWIRGQr8ClgF/D2mUflARBjSJYsbE57mR19dAf+93++QjDGmyUWTFMpVVXHPQrhfVe/HNWbX6uR3a8ORgc8ZGlpI+cqP/Q7HGGOaXDRJYYeI3ACcD7zmPVEtMbZh+WNI1703sYXXWsN4xpj4E01SOAsoAX6squtxD8q5J6ZR+WRIXjaLveYukorXwa7N/gZkjDFNLKojBVy10Yci0g8YBjwT27D80SU7hdXJffb2sDubjTFxJpqk8AGQLCJdgXeAi4GpsQzKLyJCsEv+3h6WFIwxcSaapCCqWgycCjyoqqcAg2Ibln96dO/Fd9oGgLJv7byCMSa+RJUURORg4DzgNa9fMHYh+WtIXhsWh93JZksKxph4E01S+BlwA/CS9+S03sB7sQ3LP/l52SzyrkBKKVoBpcU+R2SMMU2nzraPVPUD3HmFiu4VwNWxDMpPHbNS+Dx5LPfsSSG9x3CuCLTagyJjjPkea/uoGtJ9DA+FTua5ogGQYO2lGmPihyWFagzp6k40r9pcTFFxmc/RGGNM07GkUI19Wkz91hrHM8bEjzrPKYhIL+AqoGdkeVWdFLuw/DW4azYHBxZxVvA9+r98G1zzMQTt0RPGmNYvmj3dv4G/Aq8A4diG0zzkZiYzOG0bJ5d/DDuBzcuhw4F+h2WMMTEXTVLYo6oPxDySZibcKR8KvY718y0pGGPiQjTnFO4XkVtE5GARGVHxinlkPsvpNZRSdZej7l79uc/RGGNM04jmSGEIcAHwA/ZWH6nX3WoN6pbDMs1jkKxiz5p5pPodkDHGNIFoksIpQG9VLY11MM3JkK7ZvB3uwaDAKlI2LwZVEPE7LGOMialoqo++BNrEOpDmpm16Et+m9gUgtbwIigrrGMMYY1q+aI4UOgJLRWQ27mE7QOu+JLVCuGOVk81tuvkajzHGxFo0SeGWhk5cRCYC9+NaVf2Lqt5ZZfjlwJVACHfx52Wqurih82ts2b2GVyaFXQWfk37gD/0NyBhjYqzO6iNVfR9YCmR6ryVev1p5z3J+CDgOGAicIyIDqxT7p6oOUdVhwN3AvfWMP6YG9OjKA+Unc33ZT/iy3dF+h2OMMTFXZ1IQkTOBz4AzgDOBWSJyehTTHgMsV9UV3knqZ4GTIguo6vaIznTcVU3NxqCu2dxbfibPhn7ArKJ2fodjjDExF0310Y3AaFX9DkBEcoG3gWl1jNcVWBPRXQgcVLWQiFwJXAskUcNlriJyGXAZQPfu3aMIuXFkpybSKyedlZt2WRtIxpi4EM3VR4GKhODZHOV41V2/+b0jAVV9SFUPAH4N3FTdhFT1MVUdpaqjcnNzo5h14xnS1TWON7+wCNVmdSBjjDGNLpojhf+IyBvAM173WcCMKMYrBCIv18kD1tZS/lng4Sim26SGdklj4KJnGFhSQNFHK2gz7jK/QzLGmJiJ5kTzdcCjQD4wFHhMVX8dxbRnA31FpJeIJAFnA9MjC4hI34jOHwLLog28qQzOa89pwfc5PLiA4mUf1D2CMca0YLUeKXhXEL2hqkcB/6rPhFW1XEQmA2/gLkl9wnvG8+3AHFWdDkwWkaOAMmArcGFDFiKWBuW1YVb4AI4MfkHWuk8gHIaAPYbCGNM61ZoUVDUkIsUikq2q9T7TqqozqFLVpKpTIj7/rL7TbGoZyQksSB/LkXu+IKNsE3w7B7qN8TssY4yJiaiazgYWiMhbwK6Knqp6dcyiamb29D6W8KJHCIhSsuBlki0pGGNaqWjqQV4DbgY+AOZGvOLGuBFDmKvu9EfZwpdd43jGGNMK1XikICLvqOqRwMAoTyy3Wgf1asf9wYMZrV+TUbwGNiyCToP9DssYYxpdbUcKnUVkPDBJRIZHPmAnHh6yEykhGCDU//jK7pIF//YxGmOMiZ3azilMAa7H3V9QtU2iVv+QnarGjhjJwsU9GRwoYOeSt0g+utr77IwxpkWrMSmo6jRgmojcrKq/acKYmqWDD2jPNYGz2VkaJjX7iOZ3l50xxjSCaG5ei/uEAJAYDJA66HhmhofxzrIidpaU+x2SMcY0OrsLqx6Oz+8MQGl5mHeWbPA5GmOMaXyWFOrh0ANyyExxNW6fz/0Mirf4HJExxjSuqJKCiARFpIuIdK94xTqw5igpIcCp/RJ4K+k6bltzESXznvc7JGOMaVTRPGTnKmAD8BbuRrbXgFdjHFezNW7oINJlNwA7vnjJ52iMMaZxRXOk8DOgv6oO8h6dOURV82MdWHN1WL9c3vWeFdR2o1UhGWNal2iSwhrAHjvmSUkMsrXHsQAECVO6+DWfIzLGmMYTTYN4K4CZIvIaUFLRU1Wr3tAWN/qNOYZNq7LIke1sm/siHUZd4HdIxhjTKKI5UliNO5+QBGRGvOLW+P6dmMlIANqu/whKdvgckTHGNI46jxRU9TYAEcl0nboz5lE1cymJQb7regysfY9ELaN06RskDT3d77CMMWa/RXP10WAR+QJYCCwSkbkiMij2oTVvvcf8kB2aCsCWOS/6HI0xxjSOaKqPHgOuVdUeqtoD+AXweGzDav7GD8zjfR3BRs1iya4Mv8MxxphGEc2J5nRVfa+iQ1Vnikh6DGNqEVKTgszscx1XL95J+uYk5pSHSE4I+h2WMcbsl2iOFFaIyM0i0tN73QSsjHVgLcH4of0JE2BHSTn/Xb7J73CMMWa/RZMULgFygX8BL3mfL45lUC3FEQd2IDnBrcIZC9ZeoNf+AAAeKUlEQVT7HI0xxuy/aK4+2gpc3QSxtDgZyQmM75vDuqWfcuCi5ykbcyWJPcb4HZYxxjRYbc9ovk9Vfy4ir+CetLYPVZ0U08haiJMHpHHMiikkEKbwo0zyLCkYY1qw2o4UnvLef98UgbRU4/L7MevVQRwqC8hc+R8IhyFgLZIbY1qmGvdeqjrX+zhMVd+PfAHDmia85i8zJZGVuUcAkF2+kfLCOT5HZIwxDRfNX9oLq+l3USPH0aK1HXEqYRUA1n/6gs/RGGNMw9WYFETkHO98Qi8RmR7xeg/Y3HQhNn+HjRjMF9oXgJTlM0C/dwrGGGNahNrOKXwMrANygD9E9N8BzI9lUC1Ndmoiy9pPYOTWr8kpLaR83UISugzxOyxjjKm32s4prFLVmap6cJVzCp+ranlTBtkSZA47pfLzWqtCMsa0UNE0iDdWRGaLyE4RKRWRkIhsj2biIjJRRL4SkeUicn01w68VkcUiMl9E3hGRHg1ZiObg0NGjWBx24Sd+bQ/eMca0TNGcaP4TcA6wDEgFfgI8WNdIIhIEHgKOAwYC54jIwCrFvgBGeY/3nAbcHX3ozUubtCQ+az+Jx8uP5+ayiwiF7byCMablieqCelVdDgRVNaSqfwOOiGK0McByVV2hqqXAs8BJVab7nqoWe52fAnnRh978pBx8Gb8tP5+3d/VmToE9u9kY0/JEkxSKRSQJmCcid4vINUA0raR2xT3fuUKh168mPwZer26AiFwmInNEZM7GjRujmLU/jhnUiWDAXZr6+kJrC8kY0/JEkxQuAILAZGAX0A04LYrxpJp+1dapiMj5wCjgnuqGq+pjqjpKVUfl5uZGMWt/tEtPYmzvdgC8vmAt4bKSOsYwxpjmJZoG8VZ5H3cDt9Vj2oW4BFIhD1hbtZCIHAXcCIxX1Ra/Fz3hwGyOKriXY0tns+bda+hxrLUlaIxpOWq7ee15732Bd3XQPq8opj0b6Csivbzqp7OB6VXmMRx4FJikqt81fDGaj6PyezIxOJsusoXQwul1j2CMMc1IbUcKP/PeT2jIhFW1XEQmA2/gqp+eUNVFInI7MEdVp+OqizKAF0QEYHVLb301NyuV/2QcRudd0+m+Yy7hXVsIpLfzOyxjjIlKjUlBVdd5H08FnlfVb+s7cVWdAcyo0m9KxOej6jvNlkAGnAhzppNAmFWfvkiPIy/1OyRjjIlKNCeas4A3ReRDEblSRDrGOqiWbthhx7NFMwAoWfCyz9EYY0z06kwKqnqbqg4CrgS6AO+LyNsxj6wF69gmg3lphwDQY9unaMkOnyMyxpjo1OdpMN8B63EtpHaITTitR3l/dyommTJWz7ITzsaYliGato/+R0RmAu/gWky91GuWwtRiyLhJ7NBUAHZ9+ZLP0RhjTHTqvE8B6AH8XFXnxTqY1qRz+7Z8lDKaw0o+IHHrcjQcRuwxncaYZi6acwrXAxkicjGAiOSKSK+YR9YKbBh2FRNL7uTo4jtYvN7OKxhjmr9oqo9uAX4N3OD1SgSejmVQrcVBYw9jqXYHhNcXWFtIxpjmL5r6jFOASbh2j1DVtUBmLINqLfLapjG0WxsAXp2/lrA1p22MaeaiSQqlqqp4jdmJSDQtpBrPifmdEcIM3PouX778R7/DMcaYWkVzovl5EXkUaCMilwKXAI/HNqzW49yDutPx3Ws5Ud9j15eplBx5PslZdkWvMaZ5iuZE8+9xT0V7EegPTFHVOp+8Zpy0pATSR54BQDq7+Xra7T5HZIwxNYv2yWtvqep1qvpLVX0r1kG1NuMmns28wGAA+q1+lh0bVvockTHGVK+2prN3iMj2ml5NGWRLl5gQpPjwGwF3h3PBtJt9jsgYY6pXY1JQ1UxVzQLuA67HPUozD3d56h1NE17rcfD44/g0aSwAA797lc0FC3yOyBhjvi+a6qNjVfXPqrpDVber6sNE9zhOE0FESDnmFsIqBEVZ99JNfodkjDHfE01SCInIeSISFJGAiJwHhGIdWGs0bNQhfJzhHiExuGgm3y78yOeIjDFmX9EkhXOBM4EN3usMr59pgI6TbqNUg+zQVN748GO/wzHGmH3UeZ+CqhYAJ8U+lPjQt/8gnuh2Gw8uz2HrqiyGr97K8O5t/Q7LGGOA+j1PwTSSo0/9MbuCrvmLO19firth3Bhj/GdJwQfd2qVx/tgeAMxauYWZX23wOSJjjHEsKfhk8g/60Cm5lF8kPE/XF35IuLzM75CMMSb6pCAiY0XkXRH5r4icHMug4kG79CTuOWA+VyX8m36h5cx79WG/QzLGmFrvaO5Upde1uCa0JwK/iWVQ8WLkadeygfYAdJ13HyV7dvkckTEm3tV2pPCIiNwsIile9zbcpahnAdbMRSNIS89k5eCrAOjIZua9+AefIzLGxLvamrk4GZgHvCoiFwA/B8JAGmDVR41k5ElXslq6AtBv2WNsL9ric0TGmHhW6zkFVX0FOBZoA/wL+EpVH1DVjU0RXDxITExi05jrAGjLDha+8FufIzLGxLPazilMEpGPgHeBhcDZwCki8oyIHNBUAcaD4cdeyPKEPgAMXfMUG9ev8TkiY0y8qu1I4Q7cUcJpwF2quk1VrwWmAPZ3thFJIEDphCkApEsJy6bd6m9Axpi4VVszF0W4o4NU4LuKnqq6zOtvGtHAw05i8UfDCBZv4ol1vei0cSe9czP8DssYE2dqO1I4BXdSuZwGNoAnIhNF5CsRWS4i11cz/HAR+VxEykXk9IbMozUJnjGV48vu5O3QcH7/5ld+h2OMiUO1XX20SVUfVNVHVLXel6CKSBB4CDgOGAicIyIDqxRbDVwE/LO+02+N+h/Qi5OGdQNgxoL1fLF6q88RGWPiTSybuRgDLFfVFapaCjxLldZWVbVAVefjLnU1wDVH9yMp6L6Wp6a/aY3lGWOaVCyTQlcg8jKaQq+fqUW3dmlMHp7AXxLv4d5Nl/H5x2/6HZIxJo7EMilINf0a9LdXRC4TkTkiMmfjxtZ/i8SPDj2AcYGFACS+ezvhkB1IGWOaRiyTQiHQLaI7D1jbkAmp6mOqOkpVR+Xm5jZKcM1Zm869WZx3JgD5oYV88vbzPkdkjIkXsUwKs4G+ItJLRJJwl7FOj+H8WpX+Z9zCTlIByPn0LnbuLvE5ImNMPIhZUlDVcmAy8AawBHheVReJyO0iMglAREaLSCHuuc+PisiiWMXT0qS16cg3fS4GoL+u4JOHLqWsPORzVMaY1k5a2tUto0aN0jlz5vgdRpMI7dnJ2j8cSreyAgBe6fg/nHD57xCp7nSNMcbUTETmquqousrZk9easWBKBu0vm86mgHvmwokbHua1Zx7yOSpjTGtmSaGZS8vtQeD8aewilSJN46kFe/jnrNV+h2WMaaVqa/vINBPteo9g7UlPctX0QuZqJ2b/ewEds5I5ckBHv0MzxrQydqTQQnQZfiw3XnwqKYkBwgqT//kFX1ozGMaYRmZJoQUZ0b0tD54zgoBAfvlCdj8xiVVrv6t7RGOMiZIlhRbm6IEdeWh8mL8n/Y6xzGftX89m8/ZdfodljGklLCm0QMcdPZFv24wG4ODQXOb8+RKKS8p8jsoY0xpYUmiJgon0+p/nKUzuC8Cxe/7D6w//inJrI8kYs58sKbRQkpJFh8unsynYAYDTtj3Bi1Pvtaa2jTH7xZJCC5bUtgupF7/EDkkH4JTV/8e///WMz1EZY1oySwotXHreYMpOf5pSEkiSEEfOv5Y333vP77CMMS2UJYVWoN2gH7Dl6PsAWKGdufmtdXzwdet/7oQxpvFZUmglOh16AcuPeJgLw1PYEM7mf56ey8Jvi/wOyxjTwlhSaEX6jD+XO88aiwjsKg1x8dTZrNls9zAYY6JnSaGVOW5IZ6acMBBQzip+loUPX0DhFksMxpjoWFJohS4+tBd/7fMxv0x8gePK32H1/RN55YPZdrmqMaZOlhRaqSPO+jmbkt0jsg+R+Yx/ZxJ/eeh3fFe02+fIjDHNmSWFViqQ2YGcaz5mbe/TAciSYi7ddBfz/3gS//lsoc/RGWOaK0sKrVlKFl1+9Fd2nvo024NtATiKWYx67TgeffQBtuwq9TlAY0xzY0khDmTkn0jWtXNZ1/VYAHJkOxetvZ1z732JtxZv8Dk6Y0xzYkkhXqS3p/NPnmP78Q+zK5DB3eVnsnRXBpf+fQ6/fOFLtu+xVlaNMfY4zvgiQtaYc9EBP2DQ16VkvrKEHXvKmTa3kB1ff8iFp57EIQO6+R2lMcZHdqQQhySzE6eO7M6b1xzOuL45dJMN3Ft6O52eOYpH//EsxaXlfodojPGJJYU41jk7lb9fMobH+s0lXUroHVjPT76+nH/dcxlzv1nvd3jGGB9YUohzIsKACx9g66E3UUYiQVHOL3uRtCeP5i/TprPDzjUYE1ekpd3lOmrUKJ0zZ47fYbRKofWL2PqPS8jZsRSAUg3yCoezvu+5jJ9wDIPz2vgcoTGmoURkrqqOqrOcJQWzj/JStvznDrLnPEgQ93jPsAoHlzxIp7xenHdQD04Y2pm0JLtGwZiWxJKC2S+hwi/Y/J//o33h28wMDeXHZddVDuuQUs5FA4SjJhxBv46ZPkZpjIlWtEnB/u6ZagXzhtPhJy/A9nUM3rSJq5cn8ezsNXy3o4Qjyj7kiiWPM3tRPx5sexI9xp3LMUN7kJIY9DtsY8x+siMFE7WyUJh3lmyg7/STOKD0q8r+WzWDV+QIduVfwMTDD6NXTrqPURpjqtMsqo9EZCJwPxAE/qKqd1YZngz8HRgJbAbOUtWC2qZpSaEZ2FrAto/+QuKXT5NevnWfQR+FBjE39xT6jT+LIwZ2taMHY5oJ35OCiASBr4GjgUJgNnCOqi6OKHMFkK+ql4vI2cApqnpWbdO1pNCMlJdStng62z54jNxNs/YZdG3p5bwi4xnSNZvRvdoxpmc7RvVoR3Zaok/BGhPfmsM5hTHAclVd4QX0LHASsDiizEnArd7nacCfRES0pdVpxauEJBLzTyc3/3TY+DVbP3yMlEXPUhYKMyN8EGUon6/exuertzLxk/P5QHNYnTEU7T6WbgeOYkzvXDpnp/q9FMaYCLFMCl2BNRHdhcBBNZVR1XIRKQLaA5siC4nIZcBlAN27d49VvGZ/5Paj7am/hxN/g3y7kD/uyOOzgi3MLthC0dpvGB5YznCWw+5P4atH2b40jTnhfrycPJiyrgeR2/8QRvXpyAG5GYiI30tjTNyKZVKo7pdd9QggmjKo6mPAY+Cqj/Y/NBMziamk9BzNcbjnRQPs+rY9G18/hvT1s0nzzkFkSTE/CM7jB+XzYNXTlBQkMmH6vZSkd2ZUj7YM6pJN93YpdG+fQY/2abRPT7JkYUwTiGVSKAQim9zMA9bWUKZQRBKAbGBLDGMyPkjvOpD0n7wAqrB5OWUrP6Loqw9JLJxF9p5CALaTyjrawa5S3ly8gU8Wr2RW8pWs0Vy+0I6sC3RkZ1oe5dk9ScrpRWbnPuTltqVHuzS6tk0lMWgtthjTGGKZFGYDfUWkF/AtcDZwbpUy04ELgU+A04F37XxCKyYCOX1JzOlLzuiLXb/t6wiv+oTyTRv5TeoQZq/cwtxVW2lTVECalNBfCumPSxzs9l7rIbxAWE9bTim5nY3Sji5tUunRPo0xqWvJbNuBlKwcMjOzaJuWRJu0RNqkJdI2LYm0pKAdcRhTi1hfkno8cB/uktQnVPW3InI7MEdVp4tICvAUMBx3hHB2xYnpmtjVR/Fhz4bl7P7oz5RvXklC0SrSiwtJ0pJ9yoRUOLDkScoq/9soi5IvIV1cuT2ayDYy2KoZFHnvOyST55NOZkd6L7LTEmmblsiBgUKyUhJIyswhOaMNSSnppCQlkJwYJDUxSEpikJTEACkJQVKTgqQkBElODJCcELAEsx9UlbBCKKyEVQmFlZAq4bB6/djbP6JMWJVQSAmFytFwOeHyckLhcjQUoiyQTHkgiXAYQqoEizchJdvRcDkaDhEOuXe88qoh9iRkszW9V+V8k3dvoH3RItCQVzaEhsOohpCwG6ecROa3n0go7GJMKN/FqI3/AnXlRcOIhkDDXrfrNyP7LIoCbd0yKJy27QnalG/yynsvQhGfw7yadByzgiMJq9IuPYkXLj+kQeu7OVx9hKrOAGZU6Tcl4vMe4IxYxmBappSOfUg57d69PVRh5wbYspLwlpXsWr+cnUWb+e0BI1i1ZRert+ymaGMh6Vv2Jo4UKaMTW+kk+95L8c9dE/hq547K7p8m/YYRgeX7lCnRBEpIYg9J7NFE/hQ6medDR1QOPz34PkcE5lEmyZQHkikPJhMKJKMSBAmiEuDbxO7MSh1HUIRgQOii6xi6Zy4SCEIgCJLg3iteQEgSWZg1DvFOtyVoKYN3fAjsPQG3bx5yHV9ljGF3IAPF7aSGFH1ASmgHEEZU3c4JRTQM6votTRnKmuQDKnfOg4pn0XPPV5U7I8IhAlSUDyGE+SbQm7eSj3I77LDSr2wpk0pfqxwe8MoGCLuXhthJKtfqtZU79SzdzhPB31WWCUa8EsSNm0CYM0tvpkA7Vy7pjKQb6COFJBAmINX/mb2p7GKeDh1d2f3bhL9yXsI71Zat8GpoLDeWXV3ZfULgE/6U9GCt42zULK5a3L+yO4cibkp5uNZxAO78bgzLdW/Lw1OS3qdPoGqtepX4igexLOTmVbQ79s86sWYuTMsgApmdILMTgR4HkwlkAmdGlintDyueoaxoHSU7NlG6YzPlOzeju7cS2L2VYMk2ksqKOKRvXzqEOrCtuIytxaW03178vdklSznJlJNFMQikhkr3GT5YVvLDYMS9GSHvFeH1XaN5ZOOQyu6JgS+5KOlPtS7mFs3g9pLeld1t2c5dKbfVvm6Ao0vuZpnmVXZfnvQI/QLf1jrOzWUX8VYoo7J7QsJMTkt4q9ZxZoTGcF/ZmMruAwJrODppZq3jbNEMdpbs3ZmFCDMkcWWt4wAkse8OMEiIJAnVUNoJeI047p1X3eeaGjJOAkpqYpBgQAgIZEpSNZfIQPne9EiYIL1z0klOzCIYEESE4m25bAyVoRJACRAWV1YliBJAJUD3Dl05LqsTARGyUmN/n48lBdN6JKXBgceTCNT20/lV1R6r/oruWE/Jjs2U7txKqHQ3odLdhMt2Ey7dDWW7Ob7beIa2G8rushB7ysIMWNaLLd/1IBgqIRjeQ0K4hMRwqfu37O0d2memMq59DqGwUh5W+hanQlHtixAQoWubvfduZGsZlNQygqdDeiI7gykEvMOIYEmw2p1UpM5ZSRyYnEnAO5JptysF9rhhYW9Jwt6OSr3P2VlZHNWhAyJuZziguAubNnXeuyPzyqoE3XgSpCQhk0v69iIYcMuXosWsWHaoKyvuKGmfdwmigSA/7TOcPWmdCIoQCAh7VpzLV2VbIRBEJOgdcbmjLfFep3Yay3HtBxEMCMEAZGxMZfWO45FAAhJMIBDcW1YCCUggyKiMjszqOKhyPQRLxrBr+yQCAa98MIFgMIFAIFh5lNc2kMCSzE57V6YqlB0HUnEUGIRA4Hs72Me+9y18UOd3e1mdJRqXtX1kTGNThXAIUAhGpKfyEthT5IZ5dc+EyyvrnQG3Q8npu3ecUDls+abuebbpAYkpe7u3r3XTlID3kojPXndiGiQk7x2nIoaK4aZVaRbnFIyJSyIQrOanlZAMGR3qN61gAuT2r7tcVVld6j9OwNqpMvY4TmOMMREsKRhjjKlkScEYY0wlSwrGGGMqWVIwxhhTyZKCMcaYSpYUjDHGVGpxN6+JyEZgVQNHz6HKA3yaCYurfiyu+muusVlc9bM/cfVQ1dy6CrW4pLA/RGRONHf0NTWLq34srvprrrFZXPXTFHFZ9ZExxphKlhSMMcZUirek8P1GCpsHi6t+LK76a66xWVz1E/O44uqcgjHGmNrF25GCMcaYWlhSMMYYU6lVJgURmSgiX4nIchG5vprhySLynDd8loj0bIKYuonIeyKyREQWicjPqikzQUSKRGSe95pS3bRiEFuBiCzw5vm9JxiJ84C3vuaLyIgmiKl/xHqYJyLbReTnVco02foSkSdE5DsRWRjRr52IvCUiy7z3tjWMe6FXZpmIXBjjmO4RkaXe9/SSiLSpYdxav/MYxXariHwb8X0dX8O4tf5+YxDXcxExFYjIvBrGjck6q2nf4Nv2paqt6gUEgW+A3kAS8CUwsEqZK4BHvM9nA881QVydgRHe50zg62rimgC86sM6KwByahl+PPA67inxY4FZPnyn63E33/iyvoDDgRHAwoh+dwPXe5+vB+6qZrx2wArvva33uW0MYzoGSPA+31VdTNF85zGK7Vbgl1F817X+fhs7rirD/wBMacp1VtO+wa/tqzUeKYwBlqvqClUtBZ4FTqpS5iTgSe/zNOBIkdg+f1BV16nq597nHcASoGss59mITgL+rs6nQBsR6dyE8z8S+EZVG3on+35T1Q+ALVV6R25HTwInVzPqscBbqrpFVbcCbwETYxWTqr6pqhVPvP8UyGuMedVXDesrGtH8fmMSl7cPOBN4prHmF2VMNe0bfNm+WmNS6Aqsiegu5Ps738oy3g+oCGjfJNEBXnXVcGBWNYMPFpEvReR1ERnURCEp8KaIzBWR6p4THs06jaWzqfmH6sf6qtBRVdeB+2ED1T1r0891dwnuCK86dX3nsTLZq9p6oobqED/X1zhgg6ouq2F4zNdZlX2DL9tXa0wK1f3jr3rdbTRlYkJEMoAXgZ+r6vYqgz/HVZEMBR4E/t0UMQGHquoI4DjgShE5vMpwP9dXEjAJeKGawX6tr/rwZd2JyI1AOfCPGorU9Z3HwsPAAcAwYB2uqqYq37Y14BxqP0qI6TqrY99Q42jV9Nuv9dUak0Ih0C2iOw9YW1MZEUkAsmnYoW69iEgi7kv/h6r+q+pwVd2uqju9zzOARBHJiXVcqrrWe/8OeAl3CB8pmnUaK8cBn6vqhqoD/FpfETZUVKN5799VU6bJ1513svEE4Dz1Kp6riuI7b3SqukFVQ6oaBh6vYZ6+bGvefuBU4LmaysRyndWwb/Bl+2qNSWE20FdEenn/Ms8GplcpMx2oOEt/OvBuTT+exuLVV/4VWKKq99ZQplPFuQ0RGYP7fjbHOK50Ecms+Iw7UbmwSrHpwI/EGQsUVRzWNoEa/735sb6qiNyOLgRerqbMG8AxItLWqy45xusXEyIyEfg1MElVi2soE813HovYIs9DnVLDPKP5/cbCUcBSVS2sbmAs11kt+wZ/tq/GPpPeHF64q2W+xl3FcKPX73bcDwUgBVcdsRz4DOjdBDEdhjusmw/M817HA5cDl3tlJgOLcFdcfAoc0gRx9fbm96U374r1FRmXAA9563MBMKqJvsc03E4+O6KfL+sLl5jWAWW4f2c/xp2HegdY5r2388qOAv4SMe4l3ra2HLg4xjEtx9UxV2xjFVfZdQFm1PadN8H6esrbfubjdnidq8bmdX/v9xvLuLz+Uyu2q4iyTbLOatk3+LJ9WTMXxhhjKrXG6iNjjDENZEnBGGNMJUsKxhhjKllSMMYYU8mSgjHGmEqWFIxpZCLSRkSu8DsOYxrCkoIxjUhEgkAbXEu89RlPRMR+j8Z3thGauCYiN3pt978tIs+IyC9FZKaIjPKG54hIgfe5p4h8KCKfe69DvP4TvPbw/4m7OetO4ACv3f17vDLXichsrzG42yKmt0RE/oxrx6mbiEwVkYXi2u2/punXiIl3CX4HYIxfRGQkrhmF4bjfwufA3FpG+Q44WlX3iEhf3N2xo7xhY4DBqrrSa+lysKoO8+ZzDNDXKyPAdK8xtdVAf9xdqFd48XRV1cHeeNU+IMeYWLKkYOLZOOAl9doIEpG62thJBP4kIsOAENAvYthnqrqyhvGO8V5feN0ZuCSxGlil7hkV4B6Q0ltEHgReA96s5/IYs98sKZh4V107L+XsrVpNieh/DbABGOoN3xMxbFct8xDgd6r66D493RFF5XiqulVEhuIenHIl7oEvl0SzEMY0FjunYOLZB8ApIpLqtYB5ote/ABjpfT49onw2sE5d088X4B4dWZ0duMcqVngDuMRrLx8R6Soi33tgitfsd0BVXwRuxj020pgmZUcKJm6p6uci8hyuVcpVwIfeoN8Dz4vIBcC7EaP8GXhRRM4A3qOGowNV3Swi/xX3cPjXVfU6ERkAfOK19L0TOB9XBRWpK/C3iKuQbtjvhTSmnqyVVGM8InIrsFNVf+93LMb4xaqPjDHGVLIjBWOMMZXsSMEYY0wlSwrGGGMqWVIwxhhTyZKCMcaYSpYUjDHGVPp/qA0d4GuH928AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dY_nonlin = 100 * (td_nonlin['Y'] - 1) \n",
    "\n",
    "plt.plot(dY[:21, 2], label='linear', linestyle='-', linewidth=2.5)\n",
    "plt.plot(dY_nonlin[:21], label='nonlinear', linestyle='--', linewidth=2.5)\n",
    "plt.title(r'Consumption response to 1% monetary policy shock')\n",
    "plt.xlabel('quarters')\n",
    "plt.ylabel('% deviation from ss')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can compute the impulse response to a version of the shock scaled down to 10% of its original size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "On iteration 0\n",
      "   max error for asset_mkt is 4.22E-06\n",
      "   max error for fisher is 2.50E-04\n",
      "   max error for wnkpc is 6.15E-08\n",
      "On iteration 1\n",
      "   max error for asset_mkt is 3.31E-06\n",
      "   max error for fisher is 1.58E-08\n",
      "   max error for wnkpc is 2.56E-06\n",
      "On iteration 2\n",
      "   max error for asset_mkt is 8.70E-08\n",
      "   max error for fisher is 2.67E-10\n",
      "   max error for wnkpc is 7.56E-09\n",
      "On iteration 3\n",
      "   max error for asset_mkt is 9.90E-10\n",
      "   max error for fisher is 1.78E-12\n",
      "   max error for wnkpc is 7.87E-11\n"
     ]
    }
   ],
   "source": [
    "td_nonlin = nonlinear.td_solve(ss, block_list, unknowns, targets,\n",
    "                               rstar=ss['r']+0.1*drstar[:,2], use_saved=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8FPX9+PHXezcXkBAghPs+RE4VwuWJWhWtireo9f56VNG21lZtlar1+/Os1qtf64lHa7VepfVCVDyqcqPcEiBI5A4ECCQkm33//vhMwhKy2QWymRzv5+Oxj8zxmZn3TGbmvTOf2c+IqmKMMcbUJOB3AMYYY+o/SxbGGGNismRhjDEmJksWxhhjYrJkYYwxJiZLFsYYY2KyZNGIiMjvRORZv+MwpiESkctE5MuI/iIR6VVHy54mIv9Ty/PcY30OVKNOFiJyoYjM8v7pa0XkfRE50u+4aoOIjBGR/Mhhqvr/VLVWd7imRkTuFJFXDmB6EZH7RaTA+zwgIhKlbEcRmSwia0RERaRHlfG/EZFNIrJARAZFDD9CRN7Z3xjr0oFuTz+parqqrvA7jvqi0SYLEbkJ+DPw/4D2QDfgL8A4P+NqDEQkye8Y6rGrgTOAQ4AhwKnANVHKhoEPgLOrjhCRjsCVQC/gKeA+b3gS8Cfgl7UdeH1k+1o9oqqN7gNkAkXAuTWUScUlkzXe589AqjduDJAP/BrYAKwFLo+Y9hRgEbAd+BG42Rt+GfBlleUo0MfrnoRLWO978f0X6OAtewuwBDgsYto84DZvWVuAF4A0oAVQjDvZFHmfTsCdwCsR058OLAQKgWlA/yrzvhn4DtgKvAakRdlWl3mxPgJsBu7xhl8BLPZi+xDo7g0Xr+wGb97fAYMitsFTwEfe9vusYjpv/OHATG+6mcDhEeOmAX/0YtkOTAHaeuPSgFeAAm99ZwLtI/aH57z/44/APUCwmvUcC5QCZd42/dYb3gmY7K17LnBVDfvVV8DVEf1XAt/E2F+TvP2kR8SwkcCrXvfBwCKv+2bgd3EcA3cC//S2yXZgPnCQtz9tAFYDJ0aUj7qO3rxeB17y5rUQyKky7ZvARmAlcGOM7Xm5t99sB1YA10TMawzu2LsFWAe8DCwATosokwxsAg6tZr0rpv+dVyYPuKjKueElL9ZVwO1AoLrjlz2P3Wa4JL0Kt29+6Q17F7ihSgzfAWdUE1tN++g0ouzbcRzLXYG3vHUqAJ6Isj4PenFn7td5dX8mqu8fbycNAUk1lLkb+AZoB2TjDvI/RuxwIa9MMi457ARae+PXAkd53a2BodX9c6rZ4SZ5O/Awb8f5BHdwXQIEcSexTyOmzfMOlK5AG29HuicixvxqThCveN0HATuAE7x1+C3uJJASMe8ZuAO9De7gvTbKtrrM2x434E5szXDfnnOB/t6w24GvvPInAbOBVrjE0R/oGLENtgNH4xL2oxXbzItjC3CxN88LvP6siANqubduzbz++7xx1wD/Bpp723IY0NIb9w7wV1ySbeet9zVR1rVyG0YM+wyX5NOAQ3EH5fFRpt8KjIzozwG2x9hfq0sWWd7/vhUwAXfi7wrMqvgfxpjnnUCJ979Iwp0gVwK/9/aHq4CV8axjxLxO8bbtvXgJEHd3YjYwEUjBXQmtAE6qYXv+FOjt7RvH4I6timNoDG5fu9/bP5rh9t3XIqYfB8yPst4V0z/sTX8M7jjo541/CfgXkAH0AL4Hrqzu+GXPY/dJ3P7W2dsGh3vzPw+YHjHNIbgT9l7/I2reR6cRfd+Oeix78/kW9+Wshff/OzJyfbz/0TO4L3TN9/u8Whsn5/r2AS4C1sUosxw4JaL/JCAvYocrJiLZ4L6NjfK6f/D+8S2rzHOPna2aHW4S8EzEuBuAxRH9g4HCiP48Ik7guIN1eUSMNSWLO4DXI8YFcN+qx0TM+2cR4x8AnoqyrS4Dfqgy7H28gyxi/juB7sBxuINwFN63tohyk4B/RPSnA+W4E+HFwIwq5b8GLos4oG6PGHcd8IHXfQUu4Q+pMn17YBfQLGLYBUQk5Wjb0Ovv6sWXETHsXmBSlOnLgYMj+vt6+4DUsC/ulSwi4pzjbevuuG+PxwPn407u/wK61LAeH0X0n4b7dh/0+jO8ZbaKtY7evKZGjBsAFHvdI6vZN24DXqhue0aJ9R3gFxH7dSkRV7m4LzTb2X1ifQP4bZR5jcElixYRw17HHQ9Bb18YEDHuGmBadcevt3364PbtYuCQapaXirsa6+v1PwT8JUps1e6jcezbUY9lYDQuse/1xdhbn+m4uwZvEseXjJo+jbXOogBoG+N+ZyfcJWWFVd6wynmoaiiifyfuxAbuHvMpwCoR+UxERu9DbOsjuour6U/fszira4ixJnusn6qGvXl1jiizLqI7cv2qs7pKf3fgUREpFJFC3AEjQGdV/QR4AvdtbL2IPC0iLaubl6oWedN2qhqzZ1WcMb+M++b0D6/C+AERSfbiTAbWRsT6V9wVRjw6AZtVdXsNMUUqAiLXtSVQpN7Ruy9U9VVVHaqqJwODcCe6ubgT0mm4q42HaphF1X1rk6qWR/SD237xrGPV7Z7mHV/dgU4V29bbvr/DJelqicjJIvKNiGz2yp8CtI0oslFVSyK2wxrcVfXZItIKOBn4Ww3rvUVVd1RZl07eMlLY+7iP9r+s0Bb3jX151RGquguXjH4mIgFcgn85ynyi7aMVou3bNR3LXYFVVc5VkfrgrsTuUtXSGtcyhsaaLL7GXTafUUOZNbgdvUI3b1hMqjpTVcfhTjjv4HYWcJeKzSvKiUiHfYg5mq5RYox18tlj/bwncrrivpHsj6rLW427ldMq4tNMVb8CUNXHVHUYMBB3Gf2biGkr10lE0nG3nyrqjiL/J+DWOWbMqlqmqnep6gDcLYJTcbf3VuNOsm0j4mypqgPjXM81QBsRyYgzpoW4WxEVDvGG7TcRaYZ7UOPXuCuV1aq6DXfPe8iBzNuzr+sYaTXudlbkfpChqqd44/fYniKSivuW+xDufn0r4D3cFw2qm8bzIvAz4Fzga1WtKbbWItKiyrqswd0CLmPv4z7Wem7CnU96Rxn/Iu5uxvHATlX9urpCNeyjsdR0LK8GutXwxXgxro7ofRHpF8eyomqUyUJVt+LuoT4pImeISHMRSfa+0TzgFXsVuF1EskWkrVc+5iN+IpIiIheJSKaqlgHbcJfw4O4dDhSRQ0UkDXcJfqCuF5EuItIG943tNW/4eiBLRDKjTPc68FMROd779vJr3Enzq1qICVwl9W0iMhBARDJF5Fyve7iIjPSWuwN3oJVHTHuKiBwpIim4Sr3pqroad9I4yHvkOUlEzsfd8vhPrGBE5FgRGSwiQdz/pAwoV9W1uMrCP4lISxEJiEhvETkmyqzWAz28b4l4cX0F3CsiaSIyBFdpHe2b7UvATSLSWUQ64bb7pBriTsPdygBI9fqruh13S2gN7hZoPxFpDxyLqx84IPuxjpFmANtE5BYRaSYiQREZJCLDvfF7bE/cN/tU3K2TkIicDJwYx3LeAYYCv8Bt41ju8o7Vo3An5X96V1WvA/8rIhki0h24iRjHvfdN/nngYRHp5K3jaC/x4SWHMK4CPNpVRdR9NI51qelYnoGrQ71PRFp4/78jqsT/Ku7cMVVEoiW8mBplsgBQ1YdxO8LtuB1zNa6isOL59HtwlYXf4Z4UmeMNi8fFQJ6IbAOuxX3jQVW/x1WKTwWW4SqXDtTfcSe7Fd7nHm9ZS3AJb4V3+b/H7SlVXerF9Tjum9FpuCdKDuhSNGL+b+MqIf/hbYcFuNsD4G69PIOrnF6Fuy0Yebvk78AfcLefhuG+laGqBbgD+9feNL8FTlXVTXGE1AF3L3sb7tvUZ+w+CVyCO0lVPFX2BtAxynz+6f0tEJE5XvcFuMrQNcDbwB9U9aMo0/8VV4k5H7dN3vWGAZU/9Doqonwx7tYVuKfhiiPG4X0bPBH3f8RLfvfhrlZuxNUP1IZ9WcdK3gn4NFyl+ErcvvYs7qkjqLI9vVtdN+JOgFuAC3FPYcVaTjHuiqQnru6mJuu8ea/BJbxrveMFXD3hDtyx9CVuX3w+1vJxT6HNx13Nbcbt+5Hnz5dwdY41JZ6a9tGoajqWI7Z/H9wXiXxcnVbVebyIOzd9IlV+zxMv2Y9bqaaOiEge8D+qOtXvWGqLiEzCVczf7ncspmERkYnAQar6sxrKjMFVqHeps8Dcci/BPTLdKH70Wx37wYsxpt7zbsNeibuqr1dEpDnu6aW/+B1LIjXa21DGmMZBRK7C3UZ+X1U/9zueSCJyEu4293rcLa1Gy25DGWOMicmuLIwxxsTUaOos2rZtqz169PA7DGOMaVBmz569SVWzY5VLaLIQkbG4tn+CwLOqel+V8am4R86G4R6VPF9V87xHuxYDS72i36jqtTUtq0ePHsyaNat2V8AYYxo5EanaakK1EpYsvB+ePIlr/CofmCkik1V1UUSxK3E/ze8jIuNxzy5XPCO8XFUPTVR8xhhj4pfIOosRQK6qrvB+CPYP9n6XxDjcT+XB/VjleO+n7MYYY+qRRCaLzuzZ+Fw+ezfYVVnGawhrK65pZoCeIjJXXEN9R1ENEbla3JvwZm3cuLF2ozfGGFMpkXUW1V0hVH1ON1qZtUA3VS0QkWHAOyIy0Gs8bXdB1aeBpwFycnLsGWBjGqGysjLy8/MpKSmJXdhElZaWRpcuXUhOTo5duBqJTBb57Nliahf2btW1oky+12piJq6pZMU1lIWqzhaRipeCWA22MU1Mfn4+GRkZ9OjRA7tLvX9UlYKCAvLz8+nZs+d+zSORt6FmAn1FpKfXuuh49m4wbDJwqdd9DvCJqqrXEmwQQER64ZplthenG9MElZSUkJWVZYniAIgIWVlZB3R1lrArC1UNicgE3Ms+gsDzqrpQRO4GZqnqZNx7kV8WkVxcS47jvcmPBu4WkRCuCd9rVXVzomI1xtRvligO3IFuw4T+zkJV38O9oyBy2MSI7hLcy0yqTvcmrjnihFNVfiwsZlcoTO/sml4UZ4wxTVeTb+7j0T/dSf4jx5H8TLR34Rhjmrr0dPdFcs2aNZxzzjk+R+OPRtPcx/7qnbKVUYHF7hXxxVugWWu/QzLG1FOdOnXijTfeSOgyQqEQSUn179Tc5K8sAp13vy5568o5NZQ0xjR1eXl5DBo0CIBJkyZx1llnMXbsWPr27ctvf/vbynJTpkxh9OjRDB06lHPPPZeiIvcyxLvvvpvhw4czaNAgrr76aipa/R4zZgy/+93vOOaYY3j00UfrfsXiUP/SVx3L6jvCvSwR2LRsBpkDjvc3IGNMVHf9eyGL1myLXXAfDejUkj+cNnCfp5s3bx5z584lNTWVfv36ccMNN9CsWTPuuecepk6dSosWLbj//vt5+OGHmThxIhMmTGDiRFdte/HFF/Of//yH0047DYDCwkI+++yzWl2v2tTkk8VBvfuyQVvRTgop/3Ge3+EYY2qwaM02pq+sPw9GHn/88WRmuteNDxgwgFWrVlFYWMiiRYs44ogjACgtLWX06NEAfPrppzzwwAPs3LmTzZs3M3DgwMpkcf75e706u15p8smiTYsUvgr2pl14NhlbFsWewBjjmwGdWtar+aamplZ2B4NBQqEQqsoJJ5zAq6++ukfZkpISrrvuOmbNmkXXrl2588479/jdQ4sWLfYv+DrS5JMFwJbM/rBlNu3LVsOuIki1R2iNqY/251ZRXRs1ahTXX389ubm59OnTh507d5Kfn0+7du0AaNu2LUVFRbzxxhsN6smqJl/BDUAHV8kdQNm+aq7PwRhjGrLs7GwmTZrEBRdcwJAhQxg1ahRLliyhVatWXHXVVQwePJgzzjiD4cOH+x3qPmk07+DOycnR/X350Vdz5nL45DEArMiZSK9Tf12LkRljDsTixYvp37+/32E0CtVtSxGZrao5saa121BA3z79ubfsApZqV45JG0MvvwMyxph6xpIFkN0yjX+1OJd120pI3yBc7ndAxhhTz1idhWdQZ/c0xIIft/ociTHG1D+WLDwDO7lnpfMKdrKtpMznaIwxpn6xZOE5tF2QR5Kf5KOU37Dx07/6HY4xxtQrVmfh6d+9I80Dc2gpxSxZPdvvcIwxpl6xKwtP+8xmLBP3HFTzzQt8jsYY0xhNmzaNU089FYDJkydz3333+RxR/OzKwiMibMw4GLYvpEPJCgiVQlKK32EZYxqp008/ndNPPz2hyygvLycYDNbKvOzKIkK4wxAAUghRvGahz9EYY+qLvLw8+vfvz1VXXcXAgQM58cQTKS4uZt68eYwaNYohQ4Zw5plnsmXLFsA1OX7LLbcwYsQIDjroIL744ou95jlp0iQmTJgAwGWXXcaNN97I4YcfTq9evfZ4Z8aDDz7I8OHDGTJkCH/4wx8qh59xxhkMGzaMgQMH8vTTT1cOT09PZ+LEiYwcOZKvv/661raBXVlEyOiZA8tc97ql0+nZ7TB/AzLG7G3u32De32su02EwnBxxi2ftd/DBbdWXPfRCOOyimItdtmwZr776Ks888wznnXceb775Jg888ACPP/44xxxzDBMnTuSuu+7iz3/+M+BeYjRjxgzee+897rrrLqZOnVrj/NeuXcuXX37JkiVLOP300znnnHOYMmUKy5YtY8aMGagqp59+Op9//jlHH300zz//PG3atKG4uJjhw4dz9tlnk5WVxY4dOxg0aBB33313zHXaF5YsIvTsdwjFH6bQTEop+cFehGRMvVT4A6z6ct+mKdkafZoeR8Y1i549e3LooYcCMGzYMJYvX05hYSHHHONeyXzppZdy7rnnVpY/66yzKsvm5eXFnP8ZZ5xBIBBgwIABrF+/HnAvUZoyZQqHHea+uBYVFbFs2TKOPvpoHnvsMd5++20AVq9ezbJly8jKyiIYDHL22WfHtU77wpJFhM5t0vlOenAI35O2ySq5jamXWnWD7jFO8B0G79mflhl9mlbd4lps1ebICwsL4ypf0XT5vsy/os0+VeW2227jmmuu2aPstGnTmDp1Kl9//TXNmzdnzJgxlc2dp6Wl1Vo9RSRLFhFEhA3pB0PR93QsXgbhcgjU/kY3xhyAwy6K67bRHjoOgcvfrdUwMjMzad26NV988QVHHXUUL7/8cuVVRm056aSTuOOOO7joootIT0/nxx9/JDk5ma1bt9K6dWuaN2/OkiVL+Oabb2p1udWxZFHF+u6nceu8dizSHrxWFqZZqiULY0z1XnzxRa699lp27txJr169eOGFF2p1/ieeeCKLFy+ufNNeeno6r7zyCmPHjuWpp55iyJAh9OvXj1GjRtXqcqtjTZRX8e53a7n+766+4q3rDmdot9YHPE9jzP6zJsprz4E0UW6PzlYxuHNmZfdCa1TQGGMASxZ76dqmGS3T3N25+ZYsjDEGsDqLvYgIl7eZz7CN79B30UYIL4aA5VRj/KSqiIjfYTRoB1rlYGfBahzcchdHB+fTMbyOXRuX+x2OMU1aWloaBQUFB3yya8pUlYKCAtLS0vZ7HnZlUY3m3YdCnutet3Q63dv39TUeY5qyLl26kJ+fz8aNG/0OpUFLS0ujS5cu+z29JYtqdDs4h9C0AEkSpihvNhz9M79DMqbJSk5OpmfPnn6H0eQl9DaUiIwVkaUikisit1YzPlVEXvPGTxeRHlXGdxORIhG5OZFxVtW9fVtW4DJwyobv6nLRxhhTLyUsWYhIEHgSOBkYAFwgIgOqFLsS2KKqfYBHgPurjH8EeD9RMUYTCAhrmh8EQPsdS8HulRpjmrhEXlmMAHJVdYWqlgL/AMZVKTMOeNHrfgM4XrxHHkTkDGAF4Etb4cVtXXPlLXU7pQWr/AjBGGPqjUQmi87A6oj+fG9YtWVUNQRsBbJEpAVwC3BXTQsQkatFZJaIzKrtyq/m3Xc3T75u6fRanbcxxjQ0iUwW1T0UXfV+TrQydwGPqGpRTQtQ1adVNUdVc7Kzs/czzOp1PngEYXXhbV954M2IGGNMQ5bIp6Hyga4R/V2ANVHK5ItIEpAJbAZGAueIyANAKyAsIiWq+kQC491Dz07tuV8vZkVZW/qkjmFgXS3YGGPqoUQmi5lAXxHpCfwIjAcurFJmMnAp8DVwDvCJul/eHFVRQETuBIrqMlEABAPCrI4XMHvVFjZttCeMjTFNW8JuQ3l1EBOAD4HFwOuqulBE7haRireUP4ero8gFbgL2erzWT4M6tQRg8dpthMrDPkdjjDH+SehXZlV9D3ivyrCJEd0lwLlVp6tS/s6EBBeHQV4LtCVlYZZv3EG/Dhl+hWKMMb6ytqFqMKhjOvcnPc17KbdROu0hv8Mxxhjf2M34GvTtkEmL4CK6yQYWr53jdzjGGOMbu7KoQVIwwOpU90vuttuX+ByNMcb4x5JFDDuy3EOz2eENlBdt8jkaY4zxhyWLGFK6RP6Se4aPkRhjjH8sWcTQvt+Iyu7C5TN9jMQYY/xjySKG3j16slbbuJ613/objDHG+MSSRQwpSQFWpbg35bXZttjnaIwxxh+WLOJQ1MZVcncsX0N4Z6HP0RhjTN2z31nEYVfvk7gjP8zCcA/+tB16Nvc7ImOMqVuWLOLQpf8oXv6kHIAF63bSs30rnyMyxpi6Zbeh4tCvQwbBgHu3xYI1W32Oxhhj6l7MZCEiR3hvrkNEfiYiD4tI98SHVn+kJQfp2y4dgAU/WrIwxjQ98VxZ/B+wU0QOAX4LrAJeSmhU9dCZmbk8m/wgf8q/EN1V4wv8jDGm0YknWYS8FxKNAx5V1UeBJtdW90GtwvwkOJcOFLAhd7bf4RhjTJ2KJ1lsF5HbgJ8B74pIEEhObFj1T9u+u3/Jvel7a/bDGNO0xJMszgd2AVeq6jqgM/BgQqOqh/r0HUChtgCgfM08n6Mxxpi6FdeVBe720xcichBwKPBqYsOqf5qlJrEiqQ8AmYWLfI7GGGPqVjzJ4nMgVUQ6Ax8DlwOTEhlUfVWY2R+AzmWr0LJin6Mxxpi6E0+yEFXdCZwFPK6qZwIDExtW/SSdDgUgiXI2rbBbUcaYpiOuZCEio4GLgHe9YcHEhVR/ZfUZXtm9wSq5jTFNSDzJ4hfAbcDbqrpQRHoBnyY2rPqp18FDKNI0AEL5c32Oxhhj6k7MtqFU9XNcvUVF/wrgxkQGVV+lp6Xw57TL+H57Cs2TR3KI3wEZY0wdsYYE99GKHuN579s1ZK9P9TsUY4ypM9aQ4D4a3DkTgI3bd7FhW4nP0RhjTN2wZLGPBnZuWdltLdAaY5qKmLehRKQncAPQI7K8qp6euLDqr4GdMrk96WVyAt8T+uJwOPgJv0MyxpiEi6fO4h3gOeDfQDix4dR/mc2SOTwllwHh5Szb1MzvcIwxpk7EkyxKVPWxhEfSgBRk9Iety+i8KxfC5RBokj87McY0IfHUWTwqIn8QkdEiMrTik/DI6rHy9u6h2eaUsGX1Yp+jMcaYxIsnWQwGrgLuA/7kfR6KZ+YiMlZElopIrojcWs34VBF5zRs/XUR6eMNHiMg87/OtiJwZ7wrVhZa9h1V2r13yjY+RGGNM3YjnNtSZQC9VLd2XGXvvvXgSOAHIB2aKyGRVjWyy9Upgi6r2EZHxwP24JtEXADmqGhKRjsC3IvJvVQ3tSwyJ0vPgHErfC5Ii5exabb/kNsY0fvFcWXwLtNqPeY8AclV1hZdo/oF7216kccCLXvcbwPEiIqq6MyIxpAG6H8tPmNaZGawMuNeQNy9Y4HM0xhiTePFcWbQHlojITNxLkIC4Hp3tDKyO6M8HRkYr411FbAWygE0iMhJ4HugOXFxfrioqbEw/mH7bV9Cp+HtQBRG/QzLGmISJJ1n8YT/nXd3Zs+oVQtQyqjodGCgi/YEXReR9Vd3jJ9MicjVwNUC3bt32M8z9E2o/GLa/RwY72bZ2GS07HVSnyzfGmLoU8zaUqn4GLAEyvM9ib1gs+UDXiP4uwJpoZUQkCcgENldZ/mJgBzComtieVtUcVc3Jzs6OI6Tak9bveP5YdhHn77qD+dta1OmyjTGmrsVMFiJyHjADOBc4D5guIufEMe+ZQF8R6SkiKcB4YHKVMpOBS73uc4BPVFW9aZK85XcH+gF5cSyzzvQ++FCeK/8p07U/89dbG1HGmMYtnttQvweGq+oGABHJBqbiKqSj8uogJgAf4l6W9Lz3Poy7gVmqOhn3y/CXRSQXd0Ux3pv8SOBWESnD/Wr8OlXdtO+rlzjZGal0aJnGum0lzP/R2ogyxjRu8SSLQEWi8BQQZwOEqvoe8F6VYRMjuktwVyxVp3sZeDmeZfhpUOdM1m0rYWn+JqvkNsY0avEkiw9E5EPgVa//fKokgKbqyDbb+Gnykxy/Yy4bl71N9kGj/A7JGGMSIp4K7t8AfwWGAIcAT6vqLYkOrCEY2bcDZwb/S0vZydqva7wrZ4wxDVqNVxber7A/VNWfAG/VTUgNx8H9+rNI+jBAc2mz+kPibAXFGGManBqvLFS1HNgpIpl1FE+DIiKs6fgTALqEfqBw9aIYUxhjTMMUT0V1CTBfRJ4TkccqPokOrKFoO3z3U8Q//Pc1HyMxxpjEiSdZvAvcAXwOzI74GGDwITmsoAsA6Sve9zkaY4xJjKh1FiLysaoeDwywCu3oggEhL/tYem18mV6lS9mxcRUtsrv7HZYxxtSqmq4sOorIMcDpInJY5IuPmvrLj6rKOOysyu6VX77uYyTGGJMYNT0NNRG4Fdem08NVxilwXKKCamiGDD+aNR+2pZNsYkved36HY4wxtS5qslDVN4A3ROQOVf1jHcbU4KQmJ/FGt9t5bZlQGO7A7LJy0pLtvdzGmMYjnh/lWaKIQ+/hY/mRbHaUlvPV8nrVjJUxxhywuNp4MrGN6ZdNSpLbnB8sWOdzNMYYU7ssWdSSFqlJHN03m+aUoAvfIVS8ze+QjDGm1sTTkGBFsx/tI8ur6g+JCqqhGt95I0+suIY0LWPp173od9zFfodkjDG1Ip6XH90ArAc+wv1A713gPwmOq0EaNvwIQriK7dL57/gcjTHG1J54rix+AfRT1YJEB9PQtc5sydfNRzK6+DN6bfkvWlaCJKf5HZYxxhyweOosVgP2Krg4lfU9BYAA9qbBAAAf7UlEQVQWFLNypjX/YYxpHOJJFiuAaSJym4jcVPFJdGANVb8jz2aXugu27fPe9jkaY4ypHfEkix9w9RUpQEbEx1Sjfbts5qceBkC3jdPQ8pDPERljzIGLWWehqncBiEiG69WihEfVwBX1PBmWzqS1bmX1d9PoethP/A7JGGMOSDxPQw0SkbnAAmChiMwWkYGJD63h6nXkuZSrALBppr1u1RjT8MXzNNTTwE2q+imAiIwBngEOT2BcDVq3rt2YmzyY4l0h/ru1E4f5HZAxxhygeJJFi4pEAaCq00SkRQJjahQ+G/E0f/5kBRTA+M076dqmud8hGWPMfovraSgRuUNEenif24GViQ6soTtxUOfK7g8XWltRxpiGLZ5kcQWQDbwFvO11X57IoBqD/h0z6OZdTViyMMY0dPE8DbUFuLEOYmlURIRxB6ezbvp7nPTjLDavepI23Qf4HZYxxuyXmt7B/WdV/aWI/Bv3Zrw9qOrpCY2sERjbHQbOeRqAuf99jTbd7/I5ImOM2T81XVm87P19qC4CaYz6D85h5dud6cmPZOR9AFiyMMY0TFHrLFR1ttd5qKp+FvkBDq2b8Bq2QED4oZ17VXmf0iVs22CtuhtjGqZ4KrgvrWbYZbUcR6OVOfSsyu6VX7zmYyTGGLP/oiYLEbnAq6/oKSKTIz6fAtZceZwG5hzDOrIASM191+dojDFm/9RUZ/EVsBZoC/wpYvh24Lt4Zi4iY4FHgSDwrKreV2V8KvASMAyXgM5X1TwROQG4D9d4YSnwG1X9JK41qmeSk4Isa3MMHTa/RZ+d31JcuJFmrbL9DssYY/ZJTXUWq1R1mqqOrlJnMUdVYzal6r2K9UngZGAAcIGIVH129Epgi6r2AR4B7veGbwJOU9XBuNtgL9OApQ0eB0CShFn2pbUVZYxpeOJpSHCUiMwUkSIRKRWRchHZFse8RwC5qrpCVUuBfwDjqpQZB7zodb8BHC8ioqpzVXWNN3whkOZdhTRIgw8/mS3qWnWXJfZGWmNMwxNP21BPAOOBfwI5wCVAnzim64x7y16FfGBktDKqGhKRrUAW7sqiwtnAXFXdVXUBInI1cDVAt27d4gjJH2mpqbyTfRmz1+5ietEIpobCpCTF82yBMcbUD3GdsVQ1FwiqarmqvgAcG8dkUt2s9qWM1xT6/cA1UeJ6WlVzVDUnO7t+1wM0O/Ln/LN8DD+UNOebFfZ8gDGmYYknWewUkRRgnog8ICK/AuJpdTYf6BrR3wVYE62MiCQBmcBmr78Lri2qS1R1eRzLq9eOPbgdyUGXG62tKGNMQxNPsrgY9zTTBGAH7uR+dhzTzQT6ikhPL9mMByZXKTOZ3b/jOAf4RFVVRFoB7wK3qep/41hWvdcyLZnDe7cF4JsFyygvK/U5ImOMiV/MZOE9FVWsqttU9S5Vvcm7LRVruhAuwXwILAZeV9WFInK3iFS0K/UckCUiucBNwK3e8Am4epE7RGSe92m3H+tXr5zTcxevJt/DlNAVLJ9uv7kwxjQcorpXG4FuhMjrqnqeiMyn+oYEhyQ6uH2Rk5Ojs2bN8juMGhUUbCL9sX6kSog52eMYev1LfodkjGniRGS2qubEKlfT01C/8P6eWjshmaystsxJPYyhpTPpvvEztDyEBON5IM0YY/xV04/y1nqdZwEh73ZU5aduwmt8dvQ+GYAsCsn79tMYpY0xpn6Ip4K7JTBFRL4QketFpH2ig2rMeh95LuXqnooqmPmWz9EYY0x84qngvktVBwLXA52Az0RkasIja6Q6de7GouSBAHRZNxWi1BkZY0x9si8/I94ArMM1+Nfgn0zyU2G3EwHooBvIXzzd52iMMSa2eNqG+rmITAM+xrVAe1V9exKqoel6xHmV3eumW8OCxpj6L55HcboDv1TVeYkOpqno0bs/3wd60638B37csJGYz6wZY4zP4qmzuBVIF5HLAUQkW0R6JjyyRu7zwfcydNdf+cWW81i9eaff4RhjTI3iuQ31B+AW4DZvUDLwSiKDagpGjxzFTtIA+Mu0mD+IN8YYX8VTwX0mcDquXSi890xkJDKopmBgp0x+0t89J/D6rHxWbtrhc0TGGBNdPMmiVF2bIAogIvG0OGvicNMJ/WjJDn4ReI3NL17kdzjGGBNVPMnidRH5K9BKRK4CpgLPJDaspmFAp5Y81uEDbkx6h2HbPyVvboN8zbgxpgmIp4L7IdwrT98E+gETVfXxRAfWVPQcdyu7NBmAkg/vtB/pGWPqpXjflPeRqv5GVW9W1Y8SHVRT0r3nQcxoeyYAB5d8y7Jv7B3dxpj6J2qyEJHtIrIt2qcug2zs+px9Bzs0FQD55I92dWGMqXdqanU2Q1VbAn/GvZSoM+7VqLcA99RNeE1Dx07dmN1xPAB9ypayeNprPkdkjDF7iuc21Emq+hdV3e69Le//iO+1qmYfDDj7draqe9Cs2Zf3ouFynyMyxpjd4kkW5SJykYgERSQgIhcBdiarZW2z2/Fdt0sA6FGex4Ipk/wNyBhjIsSTLC4EzgPWe59zvWGmlg055xYKyOT7cGf+Pr+IcNjqLowx9UPMhgRVNQ8Yl/hQTGZma14Z/gITv9hJuCDAqO/WMO7Qzn6HZYwx+/Q+C1MHzjrhGNqkuzajHvnoe8rKwz5HZIwxlizqneYpSUw4tg8AeQU7mfz1Ap8jMsYYSxb10gUju3Fkyw08k/wQx0w9jZId9rMWY4y/4k4WIjJKRD4Rkf+KyBmJDKqpS00KMmFACScE59CWQua//ZDfIRljmriafsHdocqgm3BNlY8F/pjIoAzk/PR/WBnoBsBBuc9RtHWzzxEZY5qymq4snhKRO0QkzesvxD0yez5g90USLCk5mYLhNwOQSRGL3vx/PkdkjGnKamru4wxgHvAfEbkY+CUQBpoDdhuqDgw98WK+D/YFYOCql9m6aa3PERljmqoa6yxU9d/ASUAr4C1gqao+pqob6yK4pi4QDLDzyFsBaCElLHnT7v4ZY/xRU53F6SLyJfAJsAAYD5wpIq+KSO+6CrCpO+SYs1iYPMh1r3mdTWvzfI3HGNM01XRlcQ/uquJs4H5VLVTVm4CJwP/WRXAGJBCA4+4AIE3KWPHmXT5HZIxpimpKFltxVxPjgQ0VA1V1maqOj2fmIjJWRJaKSK6I3FrN+FQRec0bP11EenjDs0TkUxEpEpEn9mWFGqOBo8cyL20Ek8tH8/u1R7J6806/QzLGNDE1JYszcZXZIfaj4UARCQJPAicDA4ALRGRAlWJXAltUtQ/wCHC/N7wEuAO4eV+X21jJBX/nxrIbWFbegcc+XuZ3OMaYJqamp6E2qerjqvqUqu7Po7IjgFxVXaGqpcA/2LtBwnHAi173G8DxIiKqukNVv8QlDQMc0j2bkwa2B+DNOfnkbijyOSJjTFOSyOY+OgOrI/rzvWHVllHVEO7WV1a8CxCRq0VklojM2rix8T+g9esT+yECAQ3x/jt/8zscY0wTkshkIdUMq/qChnjKRKWqT6tqjqrmZGdn71NwDdFB7TP41UEFfJxyMzes+S253/7X75CMMU1EIpNFPtA1or8LsCZaGRFJAjIBa9eiBmcfM4xOUgDAjg/u9DcYY0yTkchkMRPoKyI9RSQF91TV5CplJgOXet3nAJ+oqr0ergadew1kbtZPATikeAbfff6OzxEZY5qChCULrw5iAvAhsBh4XVUXisjdInK6V+w5IEtEcnENFVY+XisiecDDwGUikl/Nk1RNVo+z7qJEkwHo9sn1rPr+O58jMsY0dtJYvsjn5OTorFmz/A6jzsz+15MMm/s7AFZLR1pc9yltsjv6HJUxpqERkdmqmhOrnL38qIEaNu56ZnS9AoCuupa1T59NSbH9WM8YkxiWLBqw4Zc/xNyWxwEwsGwh05/6OeFw47hSNMbUL5YsGjAJBBlw3SssTe7P8nBHJm44moemLPU7LGNMI2TJooFLTWtBh2ve4qaMB1mlHfjLtOW8NvMHv8MyxjQyliwagcy2nXjsiuNo0yIFgN+/vYCvFlvCMMbUHksWjUT3rBY8c8kwUpICnC2f0OcfR7Py+/l+h2WMaSQsWTQiw7q34dmfBLk/+RnayRYCr57Pxo3r/A7LGNMIWLJoZI4ecyKzu14GQHf9kbVPn0txcbG/QRljGjxLFo3Q0Msf5ruWYwAYUvYdc/5yGeHysL9BGWMaNEsWjZAEgvS/7u/kJvcD4IjtH/DZC7/zOSpjTENmyaKRSk5rQbtr3madtAPg2Pz/4/N3nvE5KmNMQ2XJohFr2bYzeuHrFNEcgBFzb2PW15/6HJUxpiGyZNHIdex7GOtOeoqQBpgSzuHqD4pYvHZ/3pJrjGnKLFk0AX1Gj2PGca/xi9AENu8KcMWkmazfZq83N8bEz5JFE3H4MSdy28n9AVi7tYQrJ81gR7ElDGNMfCxZNCFXHdWLC0d2I5VSrtp4L0v/dBK5y7/3OyxjTANgyaIJERHuPn0gf2n7JuOCXzE0NI+sl45l6pvPWNPmxpgaWbJoYpKCAY6+7gkWZp8CQGsp4ifzb+azh85n3cZNPkdnjKmvLFk0QcnNWzHw+ldZMeYJtpEOwLE7P6T0ycP54tP3fY7OGFMfWbJownqNuZjg9V+R22IoAN1Yz+hpF/L+E79k205rT8oYs5sliyauRXZ3+vz6Y5YMuYVSkkiSMCM3/pMLH32fb1YU+B2eMaaesGRhIBDg4LN+R9HFU8hP7sEtZVezYGsqFzzzDfe+v5hdoXK/IzTG+MyShanUpvcwOt86i6NPu5TUpACq8NfPVvDQIw+Qu8revGdMU2bJwuxBgslcPLoH7954FIM6t2SofM+tRfeT8fzRfPCvv9sjtsY0UZYsTLX6tEvnrZ8fwe+7LyIoSnvZwti5P+fDh69gXUGh3+EZY+qYJQsTVUpSgGHXPM3Kw++lmFQATi56i22PH8UH7/+LHSVlPkdojKkroto4bivk5OTorFmz/A6j0Spas4SCly6je8niymELtRfLeoxnyNgr6dWxrY/RGWP2l4jMVtWcWOXsysLEJb3TwXT/zRcsO/g6QgQBGCgrOCnvIc589GN+9ux0Ply4jpC9vtWYRinJ7wBMAxJMpu/4ewkVTiB3yv+RteRvfFB6CFtJ58vcTXyZu4kBLXdxbb8iDj/xXNpmNPM7YmNMLbHbUGb/lZexcs06Xpq3jTdm57O9JMR1wXf4bfLr5GkH5rQ/m14/uZpD+nZHRPyO1hhTjXhvQyU0WYjIWOBRIAg8q6r3VRmfCrwEDAMKgPNVNc8bdxtwJVAO3KiqH9a0LEsW/tpZGuKdOT9y5JRT6BbO3z1cU/ks7VgY/j+MOfo4mqUEfYzSGFOV78lCRILA98AJQD4wE7hAVRdFlLkOGKKq14rIeOBMVT1fRAYArwIjgE7AVOAgVY36U2JLFvWDFm9h9SfPkTbvBdqV5e8xbhb9WdXrQvofdxEHd2pDIGBXG8b4rT4ki9HAnap6ktd/G4Cq3htR5kOvzNcikgSsA7KBWyPLRpaLtjxLFvVMOEzhwg/Z8umTdN/8JQF272cn7rqf9Wm9GN6jDaN6tWFkj9YM6NyKoCUPY+pcvMkikRXcnYHVEf35wMhoZVQ1JCJbgSxv+DdVpu2cuFBNrQsEaDX4ZFoNPpmyTSvJm/I42ctepyysLNPOaHEZUxev59PFa/g69QZm0JU1rYYR7Hkk3Q85ikHd2pEctIf1jKkvEpksqvuaWPUyJlqZeKZFRK4Grgbo1q3bvsZn6khy2570uvBhKPtf1i7+mv8t7sP0lQVMX7GZ7O2LaCeFtKMQts6HeZPYNTeZufQlv+VQpMcRdBl8NIN7diQt2eo7jPFLIpNFPtA1or8LsCZKmXzvNlQmsDnOaVHVp4Gnwd2GqrXITWIkN6PjkOO4ELhwZDdUlbW5WSz/4lwy1k2nXalrrDBVyhjBIkZsXwTzX2HXd0mMLn+avl07cVi31vTMSqNH6xS6tWtN+4w0q/swpg4kMlnMBPqKSE/gR2A8cGGVMpOBS4GvgXOAT1RVRWQy8HcReRhXwd0XmJHAWI0PRIROfQ+Fvs+6AdvXs3nxpxQumkbztd/QYddKAFZpezaH0pi+cjPTV26mi2zg85RfsYYspmt7tqR2Zmd6d7RNL5q170OrzgfRrUM7OrVKI8luZRlTKxKWLLw6iAnAh7hHZ59X1YUicjcwS1UnA88BL4tILu6KYrw37UIReR1YBISA62t6Eso0EhntaTNiPG1GjHf9OzaxdelnbF9XyCVl3Zm+YjO5G4voznoConRhE11kE5QthC24z3I36UbNZFjZI7Ru3YbuWS3ontWcw4J5NGvZhuatssnIzKJ1i1RaN08hIy3Jrk6MicF+lGcalLLyMOtXzCc07zW0YAUp21aRWbKa9PD2Pcpt1eYcsuvZyv5USlmadlllf0gDFJJOoaZTSAZFgQyKkzKZ3PoSQhldaN08mdYtUuhTvpL05mmktWxLavN0UtNakJqaRlpykLTkAM2Sg1530J7mqmPhsFKuSnlYCVf8DUN5OEyovIxwqJxwOEQ4FCJcHqIsKY2wpFROEyjaiJZuJ1weQstdWS0PEQ6Xo2E3rDg1ix0tuhFWKA8rqUX5ZGxdQri8HLQcDYchHELD5RAOQbickmA6y7J/gnrTJO/awsD1/0LD5YiGQcsh7D6V/RpmSvYVFAeaEVYlHFbGrX+ctPIdiFaUCxPwygZww/7V4lwWpwwiHFYeu+AwOrXa91YT6sPTUMbUuuRggC59D4G+h+w5ongL4YKVFP64lO1rv2dbUTHXZvVmVcEO8gp2klywZI/iSRKmLdtoK9t2DwzBffk/ZZUmVw76OOXX9A6s3WPakAYoIYVi73N72UW8Hx5JSlKAtKQAlwQ/ZChLKA+kUhZsRjiYigaSQQJoIAkkwI9p/VjU8nCSAkIwIHQqzaPPznloIAkJBCEQRAJJEEiq/PV7WbA5K9scWRlHSmgHPQu/2iM28Z4NifzBfG7rIwkF0lBVFDho40cEy4sRFNVyRBU0DBV/CZObMYpNqV3dSVKV/oWf0a54BaLlEA4juJOWaNgN03JyUwcyo/kYd+JW6L9zDkfumOLKeCe3ivIBLSdAmI2SxYNpEypPrO1C67i77EECGiaIKxMgTFAjuiln7K772ELLynX8MvVG2rGFIGGCUv0X4AmlN/Cf8OjK/ieSH+XU4PRqy1Z4JXQ8t4eurOz/WfAj7kl+ocZpcsOduKp0d5VrH8nn56lP1DgNwK/yx7CJzMr+21M/2nP/rMbz20YwI9wJgJ2lib35YsnCNA7NWhPo0po2XYbSxhs0OGK0Fg+maFlHirduZNf2TYSKNlG+YzNSvJlgyRZSSgtpVlZIn+5dSStOY8vOUgp3ltFatu+1qCQJk04J6ZQAkIJrqr00FKY0FKZv8mKOC37l2h6I0or7q6Fj+XeoR2X/BcHPuCn5uRpXcUW4A78qbV/Z31PW8mnq72JtGYaXPMlGWlf2z0h9kHZS8ztJbiidwL/Dh1f2P578b04KflPDFLA+tJWPQwMq+7sGl3JU8sc1TrM83JEV23dU9ifLdgakLt+7YJWLtiT2PDGmECJFaj5ZBtizkctwHO2oJnKaUGX6c592GcmkBNwDGwERNpe0pVxTCBNAJYgihMWly7AEUAJ0bN2eUc3aEBAhLTmx9XOWLEyTIM1akT7kVNJjlIs8XasqJd+/zKYt6yndXkBo1w7CpcWES4vRsp1QVoyWlTCm40h6Nu9LSVmYkrJy2qzqwvrtXUkK7yI5vIsU3UVQQ963YnciSW+WSs/UFoTCYcJhaFUWgAR9MYz8QaRINc+gV6NFSoBMSSYYEAICyeXJlROGkYiTnHhrFaB5ixYMbNGSYEAQEbJ2tWXDjvaV1wNh2X2SU3H9W5Lb89NOHQmIEBRoXZ7C4rUjwRuPBNFAEK3SfVmPfpSlZBIUIRAQ8n64kDXhYpAgBAIQSPauzoLesCDjOhzFCS17e+skZBVMYGHxWUggiUAgCMEkAhJEgm46CSQxPKMj77Y5iGBACIqQtGsQa4rOJRAMIsEkAoEkAsEggWCQYDAJCSbRPimFJRntEYGgiLvSCV1YGUdFjFVPvu/t9V+YG/P/9Ps4/pe1xeosjKlLkbd8ghGni7Ji2LXdu5cdirivHXJlAYLJkNV79zShXejmFXvMGvZMBhoOE2jbF0lK3t2YY+FqV0oC7oPs7paAyygpLSApdfeMykPubyC45z0u0+BZnYUx9ZGI+2ZZVXIz99kXSalIu/67Zx3vdK26xi5TVdBOFU2dPYRujDEmJksWxhhjYrJkYYwxJiZLFsYYY2KyZGGMMSYmSxbGGGNismRhjDEmpkbzozwR2QisOoBZtAU21VI4tcni2jcW176xuPZNY4yru6pmxyrUaJLFgRKRWfH8irGuWVz7xuLaNxbXvmnKcdltKGOMMTFZsjDGGBOTJYvdnvY7gCgsrn1jce0bi2vfNNm4rM7CGGNMTHZlYYwxJiZLFsYYY2JqUslCRMaKyFIRyRWRW6sZnyoir3njp4tIjzqIqauIfCoii0VkoYj8opoyY0Rkq4jM8z4TEx1XxLLzRGS+t9y93i4lzmPeNvtORIYmOJ5+EdthnohsE5FfVilTZ9tLRJ4XkQ0isiBiWBsR+UhElnl/W0eZ9lKvzDIRubQO4npQRJZ4/6e3RaRVlGlr/J8nIK47ReTHiP/XKVGmrfH4TUBcr0XElCci86JMm8jtVe35wZd9TFWbxAcIAsuBXkAK8C0woEqZ64CnvO7xwGt1EFdHYKjXnQF8X01cY4D/+LTd8oC2NYw/BXgf9+6dUcD0Ov6frsP9qMiX7QUcDQwFFkQMewC41eu+Fbi/munaACu8v6297tYJjutEIMnrvr+6uOL5nycgrjuBm+P4X9d4/NZ2XFXG/wmY6MP2qvb84Mc+1pSuLEYAuaq6QlVLgX8A46qUGQe86HW/ARwvkth3SKrqWlWd43VvBxYDnRO5zFo2DnhJnW+AViLSsY6WfTywXFUP5Jf7B0RVPwc2VxkcuR+9CJxRzaQnAR+p6mZV3QJ8BIxNZFyqOkVVvfej8g3QpbaWdyBxxSme4zchcXnngPOAV2trefGq4fxQ5/tYU0oWnYHVEf357H1SrizjHVRbgaw6iQ7wbnsdBkyvZvRoEflWRN4XkYF1FRPulc5TRGS2iFxdzfh4tmuijCf6AezX9gJor6prwR3sQLtqyvi53QCuwF0RVifW/zwRJni3x56PckvFz+11FLBeVZdFGV8n26vK+aHO97GmlCyqu0Ko+txwPGUSQkTSgTeBX6rqtiqj5+ButRwCPA68UxcxeY5Q1aHAycD1InJ0lfG+bDMRSQFOB/5ZzWg/t1e8/NzXfg+EgL9FKRLrf17b/g/oDRwKrMXd8qnKt+0FXEDNVxUJ314xzg9RJ6tm2H5vs6aULPKByDfVdwHWRCsjIklAJvt3ybxPRCQZtyP8TVXfqjpeVbepapHX/R6QLCJtEx2Xt7w13t8NwNu42wGR4tmuiXAyMEdV11cd4ef28qyvuBXn/d1QTRlftptXyXkqcJF6N7ariuN/XqtUdb2qlqtqGHgmyvL82l5JwFnAa9HKJHp7RTk/1Pk+1pSSxUygr4j09L6VjgcmVykzGah4YuAc4JNoB1Rt8e6HPgcsVtWHo5TpUFF3IiIjcP+3gkTG5S2rhYhkVHTjKkgXVCk2GbhEnFHA1orL4wSL+m3Pr+0VIXI/uhT4VzVlPgROFJHW3m2XE71hCSMiY4FbgNNVdWeUMvH8z2s7rsg6rjOjLC+e4zcRfgIsUdX86kYmenvVcH6o+30sETX49fWDe3Lne9xTFb/3ht2NO3gA0nC3NXKBGUCvOojpSNyl4XfAPO9zCnAtcK1XZgKwEPcEyDfA4XW0vXp5y/zWW37FNouMTYAnvW06H8ipg7ia407+mRHDfNleuIS1FijDfZO7ElfP9TGwzPvbxiubAzwbMe0V3r6WC1xeB3Hl4u5hV+xnFU/+dQLeq+l/nuC4Xvb2ne9wJ8GOVePy+vc6fhMZlzd8UsV+FVG2LrdXtPNDne9j1tyHMcaYmJrSbShjjDH7yZKFMcaYmCxZGGOMicmShTHGmJgsWRhjjInJkoUxdUREWonIdX7HYcz+sGRhTB0QkSDQCtey8b5MJyJix6nxne2ExlRDRH7vvTthqoi8KiI3i8g0EcnxxrcVkTyvu4eIfCEic7zP4d7wMd67CP6O+9HZfUBv770HD3plfiMiM71G9O6KmN9iEfkLrp2rriIySUQWiHtvwq/qfouYpi7J7wCMqW9EZBiuOYnDcMfIHGB2DZNsAE5Q1RIR6Yv7NXCON24EMEhVV3qthg5S1UO95ZwI9PXKCDDZa4TuB6Af7he313nxdFbVQd501b60yJhEsmRhzN6OAt5Wr/0kEYnVBlEy8ISIHAqUAwdFjJuhqiujTHei95nr9afjkscPwCp17wcB99KaXiLyOPAuMGUf18eYA2bJwpjqVdcOTojdt27TIob/ClgPHOKNL4kYt6OGZQhwr6r+dY+B7gqkcjpV3SIih+BeZnM97kU8V8SzEsbUFquzMGZvnwNnikgzr0XR07zhecAwr/uciPKZwFp1TWxfjHsFaHW2416NWeFD4ArvXQWISGcR2eslNl7z6gFVfRO4A/f6T2PqlF1ZGFOFqs4RkddwLXyuAr7wRj0EvC4iFwOfREzyF+BNETkX+JQoVxOqWiAi/xWRBcD7qvobEekPfO21qF4E/Ax3KytSZ+CFiKeibjvglTRmH1mrs8bEICJ3AkWq+pDfsRjjF7sNZYwxJia7sjDGGBOTXVkYY4yJyZKFMcaYmCxZGGOMicmSxf9vrw4EAAAAAAT5W4+wQEkEwJIFACtnVwVwAWuwowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dY_nonlin = 100 * (td_nonlin['Y'] - 1) \n",
    "\n",
    "plt.plot(0.1*dY[:21, 2], label='linear', linestyle='-', linewidth=2.5)\n",
    "plt.plot(dY_nonlin[:21], label='nonlinear', linestyle='--', linewidth=2.5)\n",
    "plt.title(r'Consumption response to 0.1% monetary policy shock')\n",
    "plt.xlabel('quarters')\n",
    "plt.ylabel('% deviation from ss')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the shock is this small, the linear and nonlinear impulse responses agree exactly. This is a good check, amid a highly complex model, that we didn't make any mistakes."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
